Neonate Data Analysis Colorado 2022

Robin Bedard & Savannah Weaver

Packages

library("tidyverse") # for dplyr workflow + ggplot
library("rmdformats") # pretty Rmd format
library("lme4") # LMM
library("lmerTest") # LMM stats
library("emmeans") # pairwise comparisons
library("broom") # clean stats presentation
library("AICcmodavg") # model selection
library("MuMIn") # GLMM R-squared
library("effects") # partial regression plots
library("ggpubr") # figure arranging & exporting

Fig Theme

my_ggplot_theme <- theme_classic() + # start with a preset theme to make things easy
                     theme(axis.title = element_text(size = 10), # text size for axis labels
                           axis.text = element_text(size = 8, # text size for axis values
                                                    color = "black", # need to specifically set color
                                                    family = "sans"), # and family
                           legend.title = element_text(size = 8), # text size for the legend title
                           legend.text = element_text(size = 8), # text size for legend items
                           text = element_text(color = "black", # color for all text (except axis values)
                                               family = "sans") # family
                           )

Now, apply the theme. This sets the theme for the entire Rmd document.

theme_set(my_ggplot_theme)

Acronyms

DAB = day after birth PS = post-shed / 1 week after birth

Load & Tidy Data

DAB osmolality

This is the osmolality of each individual neonate the day after birth.

neonate_dab_osmol <- read_rds("data/neonate_DAB_osmolality.RDS")
head(neonate_dab_osmol)
## # A tibble: 6 × 3
## # Groups:   CEWL_Blood_Collect_Date [1]
##   CEWL_Blood_Collect_Date Indiv_ID Plasma_Osmol_Rep_mean
##   <date>                  <fct>                    <dbl>
## 1 2022-08-28              300                        306
## 2 2022-08-28              301                        312
## 3 2022-08-28              302                        304
## 4 2022-08-28              303                        313
## 5 2022-08-28              304                        293
## 6 2022-08-28              305                        278

DAB body size

This is the body size and body temperature at time of CEWL measurements for each individual neonate.

neonate_dab_all0 <- read_rds("data/neonate_DAB_body_size.RDS")

neonate_dab_all <- neonate_dab_all0 %>% 
  # fix a typo of where the decimal was put
  mutate(SVL_cm = case_when(SVL_cm == 2.4 ~ SVL_cm*10, TRUE ~ SVL_cm)) %>% 
  rename(Date_Born_Shed = Date_Born) %>% 
  mutate(timept = "DAB")

# robin removed this one...
#neonate_dab_all[77, ]

summary(neonate_dab_all)
##     Indiv_ID  Date_Born_Shed       Tail_Length_cm      SVL_cm     
##  300    : 1   Min.   :2022-08-26   Min.   :1.200   Min.   :20.50  
##  301    : 1   1st Qu.:2022-08-30   1st Qu.:1.600   1st Qu.:22.00  
##  302    : 1   Median :2022-09-01   Median :1.900   Median :22.55  
##  303    : 1   Mean   :2022-08-31   Mean   :1.829   Mean   :22.67  
##  304    : 1   3rd Qu.:2022-09-04   3rd Qu.:2.000   3rd Qu.:23.15  
##  305    : 1   Max.   :2022-09-05   Max.   :2.500   Max.   :25.70  
##  (Other):76                                                       
##  Tail_SVL_Ratio        Mass_g     Sex      Mother_ID  Treatment   Tb_CEWL_c    
##  Min.   :0.05330   Min.   :11.6   F:27   101    : 5   C:36      Min.   :25.80  
##  1st Qu.:0.06825   1st Qu.:13.5   M:55   124    : 5   W:46      1st Qu.:27.23  
##  Median :0.08480   Median :14.9          125    : 5             Median :28.25  
##  Mean   :0.09035   Mean   :15.4          126    : 5             Mean   :28.23  
##  3rd Qu.:0.09150   3rd Qu.:16.5          127    : 5             3rd Qu.:29.00  
##  Max.   :0.87500   Max.   :23.0          128    : 5             Max.   :31.30  
##                                          (Other):52                            
##  Pee    CEWL_Blood_Collect_Date    timept         
##  N:39   Min.   :2022-08-28      Length:82         
##  Y:43   1st Qu.:2022-08-31      Class :character  
##         Median :2022-09-02      Mode  :character  
##         Mean   :2022-09-01                        
##         3rd Qu.:2022-09-05                        
##         Max.   :2022-09-06                        
## 

PS lumped

neonate_postshed <- read_rds("data/neonate_PS_all_data.RDS")

PS osmolality

first, separate out and get mean osmolality per individual

postshed_osmol <- neonate_postshed %>% 
  select(CEWL_Blood_Collect_Date, Indiv_ID, 
         Plasma_Osmol_Rep1, Plasma_Osmol_Rep2, 
         Plasma_Osmol_Rep3, Plasma_Osmol_Rep4) %>% 
  pivot_longer(3:6, names_to = "rep", values_to = "osml") %>% 
  filter(complete.cases(osml)) %>% 
  group_by(CEWL_Blood_Collect_Date, Indiv_ID) %>% 
  summarise(Plasma_Osmol_Rep_mean = mean(osml))

head(postshed_osmol)
## # A tibble: 6 × 3
## # Groups:   CEWL_Blood_Collect_Date [1]
##   CEWL_Blood_Collect_Date Indiv_ID Plasma_Osmol_Rep_mean
##   <date>                  <fct>                    <dbl>
## 1 2022-09-03              303                       352 
## 2 2022-09-03              304                       345 
## 3 2022-09-03              344                       349 
## 4 2022-09-03              345                       332.
## 5 2022-09-03              346                       344 
## 6 2022-09-03              347                       428

PS body size

neonate_postshed_clean <- neonate_postshed %>% 
  select(Indiv_ID, Date_Born_Shed = Date_Shed,
         CEWL_Blood_Collect_Date, 
         Tail_Length_cm, SVL_cm, Tail_SVL_Ratio, 
         Mass_g, Sex, Pee, Tb_CEWL_c, 
         Mother_ID, Treatment) %>% 
  mutate(timept = "PS") %>% 
  # fix a typo
  mutate(Treatment = case_when(Mother_ID == 128 ~ "W", TRUE ~ Treatment)) %>% 
  mutate(Treatment = factor(Treatment))
summary(neonate_postshed_clean)
##     Indiv_ID  Date_Born_Shed       CEWL_Blood_Collect_Date Tail_Length_cm
##  303    : 1   Min.   :2022-09-02   Min.   :2022-09-03      Min.   :1.50  
##  304    : 1   1st Qu.:2022-09-02   1st Qu.:2022-09-03      1st Qu.:1.65  
##  320    : 1   Median :2022-09-03   Median :2022-09-04      Median :1.90  
##  344    : 1   Mean   :2022-09-03   Mean   :2022-09-04      Mean   :1.88  
##  345    : 1   3rd Qu.:2022-09-04   3rd Qu.:2022-09-05      3rd Qu.:2.05  
##  346    : 1   Max.   :2022-09-04   Max.   :2022-09-05      Max.   :2.30  
##  (Other):29                                                              
##      SVL_cm      Tail_SVL_Ratio        Mass_g      Sex    Pee   
##  Min.   :22.40   Min.   :0.06220   Min.   :11.10   F:13   N:29  
##  1st Qu.:23.15   1st Qu.:0.06855   1st Qu.:12.05   M:22   Y: 6  
##  Median :23.80   Median :0.08160   Median :13.00                
##  Mean   :23.87   Mean   :0.07881   Mean   :13.42                
##  3rd Qu.:24.00   3rd Qu.:0.08665   3rd Qu.:14.40                
##  Max.   :26.50   Max.   :0.09750   Max.   :19.00                
##                                                                 
##    Tb_CEWL_c       Mother_ID Treatment    timept         
##  Min.   :26.50   125    :5   C:19      Length:35         
##  1st Qu.:27.75   128    :5   W:16      Class :character  
##  Median :28.50   133    :5             Mode  :character  
##  Mean   :28.63   134    :5                               
##  3rd Qu.:29.55   135    :5                               
##  Max.   :30.60   129    :4                               
##                  (Other):6

CEWL

all_data_CEWL0 <- read_rds("data/CEWL_data_all.RDS")

all_data_CEWL <- all_data_CEWL0 %>% 
  rename(CEWL_Blood_Collect_Date = date)

summary(all_data_CEWL)
##  CEWL_Blood_Collect_Date    Indiv_ID     CEWL_g_m2h      msmt_temp_C   
##  Min.   :2022-08-22      103    :  4   Min.   : 4.043   Min.   :22.86  
##  1st Qu.:2022-08-29      115    :  4   1st Qu.:12.500   1st Qu.:24.95  
##  Median :2022-09-02      101    :  3   Median :15.350   Median :25.68  
##  Mean   :2022-08-31      102    :  3   Mean   :16.404   Mean   :25.49  
##  3rd Qu.:2022-09-05      104    :  3   3rd Qu.:20.636   3rd Qu.:26.12  
##  Max.   :2022-09-06      105    :  3   Max.   :30.150   Max.   :27.40  
##                          (Other):199                                   
##  msmt_RH_percent      reps          e_s_kPa         e_a_kPa     
##  Min.   :32.92   Min.   :3.000   Min.   :2.784   Min.   :1.141  
##  1st Qu.:39.16   1st Qu.:4.000   1st Qu.:3.157   1st Qu.:1.265  
##  Median :41.27   Median :4.000   Median :3.296   Median :1.342  
##  Mean   :41.14   Mean   :4.174   Mean   :3.265   Mean   :1.341  
##  3rd Qu.:42.76   3rd Qu.:4.000   3rd Qu.:3.385   3rd Qu.:1.409  
##  Max.   :49.68   Max.   :8.000   Max.   :3.648   Max.   :1.693  
##                                                                 
##     VPD_kPa     
##  Min.   :1.488  
##  1st Qu.:1.836  
##  Median :1.927  
##  Mean   :1.923  
##  3rd Qu.:2.025  
##  Max.   :2.324  
## 

Mother Data!

mother_values <- read_rds("data/mother_values.RDS")
neo_counts <- read_csv("data/neonate_count.csv", 
                       col_types = "fddd") %>% 
  mutate(
    total_embryos = (live_neonates+dead_at_birth+slugs),
    prop_live = live_neonates/total_embryos
  )

summary(neo_counts)
##    Mother_ID  live_neonates    dead_at_birth        slugs       
##  133    : 1   Min.   : 5.000   Min.   :0.0000   Min.   :0.0000  
##  112    : 1   1st Qu.: 6.500   1st Qu.:0.0000   1st Qu.:0.0000  
##  115    : 1   Median : 8.000   Median :0.0000   Median :0.0000  
##  119    : 1   Mean   : 8.105   Mean   :0.4737   Mean   :0.3684  
##  126    : 1   3rd Qu.: 9.500   3rd Qu.:1.0000   3rd Qu.:1.0000  
##  114    : 1   Max.   :13.000   Max.   :2.0000   Max.   :2.0000  
##  (Other):13                                                     
##  total_embryos      prop_live     
##  Min.   : 6.000   Min.   :0.6250  
##  1st Qu.: 7.500   1st Qu.:0.8229  
##  Median : 9.000   Median :1.0000  
##  Mean   : 8.947   Mean   :0.9109  
##  3rd Qu.:10.000   3rd Qu.:1.0000  
##  Max.   :16.000   Max.   :1.0000  
## 
neo_counts %>% 
  filter(slugs > 0 | dead_at_birth > 0)
## # A tibble: 9 × 6
##   Mother_ID live_neonates dead_at_birth slugs total_embryos prop_live
##   <fct>             <dbl>         <dbl> <dbl>         <dbl>     <dbl>
## 1 133                  13             1     2            16     0.812
## 2 126                  10             1     0            11     0.909
## 3 125                   8             2     0            10     0.8  
## 4 128                   8             2     0            10     0.8  
## 5 129                   8             0     1             9     0.889
## 6 134                   8             0     1             9     0.889
## 7 105                   6             1     1             8     0.75 
## 8 117                   5             2     1             8     0.625
## 9 131                   5             0     1             6     0.833
neo_counts %>% 
  filter(slugs > 0 | dead_at_birth > 0) %>% 
  summary()
##    Mother_ID live_neonates    dead_at_birth     slugs        total_embryos   
##  133    :1   Min.   : 5.000   Min.   :0     Min.   :0.0000   Min.   : 6.000  
##  126    :1   1st Qu.: 6.000   1st Qu.:0     1st Qu.:0.0000   1st Qu.: 8.000  
##  125    :1   Median : 8.000   Median :1     Median :1.0000   Median : 9.000  
##  128    :1   Mean   : 7.889   Mean   :1     Mean   :0.7778   Mean   : 9.667  
##  129    :1   3rd Qu.: 8.000   3rd Qu.:2     3rd Qu.:1.0000   3rd Qu.:10.000  
##  134    :1   Max.   :13.000   Max.   :2     Max.   :2.0000   Max.   :16.000  
##  (Other):3                                                                   
##    prop_live     
##  Min.   :0.6250  
##  1st Qu.:0.8000  
##  Median :0.8125  
##  Mean   :0.8120  
##  3rd Qu.:0.8889  
##  Max.   :0.9091  
## 

Join/Arrange Dataframes

DAB

All DAB: missing CEWL for neonate 332

dab_all <- neonate_dab_all %>% 
  left_join(neonate_dab_osmol, 
            by = c('Indiv_ID', 'CEWL_Blood_Collect_Date')) %>% 
  left_join(all_data_CEWL, 
            by = c('Indiv_ID', 'CEWL_Blood_Collect_Date')) %>% 
  # remove one erroneous point
  filter(Plasma_Osmol_Rep_mean < 400) %>% 
  left_join(mother_values, by = "Mother_ID") %>% 
  left_join(neo_counts, by = "Mother_ID") 

summary(dab_all)
##     Indiv_ID  Date_Born_Shed       Tail_Length_cm      SVL_cm     
##  300    : 1   Min.   :2022-08-26   Min.   :1.200   Min.   :20.50  
##  301    : 1   1st Qu.:2022-08-30   1st Qu.:1.600   1st Qu.:22.00  
##  302    : 1   Median :2022-09-01   Median :1.900   Median :22.50  
##  303    : 1   Mean   :2022-08-31   Mean   :1.835   Mean   :22.66  
##  304    : 1   3rd Qu.:2022-09-04   3rd Qu.:2.000   3rd Qu.:23.20  
##  305    : 1   Max.   :2022-09-05   Max.   :2.500   Max.   :25.70  
##  (Other):75                                                       
##  Tail_SVL_Ratio        Mass_g      Sex      Mother_ID  Treatment
##  Min.   :0.05330   Min.   :11.60   F:26   101    : 5   C:35     
##  1st Qu.:0.06870   1st Qu.:13.50   M:55   124    : 5   W:46     
##  Median :0.08520   Median :14.90          125    : 5            
##  Mean   :0.09071   Mean   :15.41          126    : 5            
##  3rd Qu.:0.09170   3rd Qu.:16.50          127    : 5            
##  Max.   :0.87500   Max.   :23.00          128    : 5            
##                                           (Other):51            
##    Tb_CEWL_c     Pee    CEWL_Blood_Collect_Date    timept         
##  Min.   :25.80   N:39   Min.   :2022-08-28      Length:81         
##  1st Qu.:27.30   Y:42   1st Qu.:2022-08-31      Class :character  
##  Median :28.30          Median :2022-09-02      Mode  :character  
##  Mean   :28.25          Mean   :2022-09-01                        
##  3rd Qu.:29.00          3rd Qu.:2022-09-05                        
##  Max.   :31.30          Max.   :2022-09-06                        
##                                                                   
##  Plasma_Osmol_Rep_mean   CEWL_g_m2h      msmt_temp_C    msmt_RH_percent
##  Min.   :263.0         Min.   : 4.043   Min.   :23.56   Min.   :32.92  
##  1st Qu.:283.7         1st Qu.: 9.841   1st Qu.:25.28   1st Qu.:38.85  
##  Median :293.0         Median :12.359   Median :26.00   Median :40.84  
##  Mean   :295.7         Mean   :11.995   Mean   :25.78   Mean   :40.82  
##  3rd Qu.:304.0         3rd Qu.:13.926   3rd Qu.:26.48   3rd Qu.:42.57  
##  Max.   :355.0         Max.   :23.808   Max.   :27.20   Max.   :49.42  
##                        NA's   :1        NA's   :1       NA's   :1      
##       reps          e_s_kPa         e_a_kPa         VPD_kPa     
##  Min.   :3.000   Min.   :2.904   Min.   :1.141   Min.   :1.566  
##  1st Qu.:4.000   1st Qu.:3.219   1st Qu.:1.294   1st Qu.:1.883  
##  Median :4.000   Median :3.360   Median :1.353   Median :1.959  
##  Mean   :4.188   Mean   :3.320   Mean   :1.354   Mean   :1.966  
##  3rd Qu.:4.000   3rd Qu.:3.455   3rd Qu.:1.416   3rd Qu.:2.077  
##  Max.   :8.000   Max.   :3.606   Max.   :1.693   Max.   :2.324  
##  NA's   :1       NA's   :1       NA's   :1       NA's   :1      
##   Mother_CEWL     Mother_mass      Mother_SVL     Mother_osml   
##  Min.   :10.27   Min.   :181.7   Min.   :72.00   Min.   :281.7  
##  1st Qu.:16.96   1st Qu.:281.2   1st Qu.:75.50   1st Qu.:304.0  
##  Median :19.90   Median :319.7   Median :80.50   Median :311.7  
##  Mean   :19.25   Mean   :337.7   Mean   :81.24   Mean   :312.7  
##  3rd Qu.:21.25   3rd Qu.:379.0   3rd Qu.:86.50   3rd Qu.:320.8  
##  Max.   :27.82   Max.   :570.9   Max.   :91.00   Max.   :357.2  
##                                                                 
##  Mother_Days_in_Treatment live_neonates    dead_at_birth        slugs       
##  Min.   : 2.000           Min.   : 5.000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.: 5.000           1st Qu.: 7.000   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median : 7.000           Median : 8.000   Median :0.0000   Median :0.0000  
##  Mean   : 8.333           Mean   : 8.284   Mean   :0.4198   Mean   :0.3457  
##  3rd Qu.:14.000           3rd Qu.:10.000   3rd Qu.:1.0000   3rd Qu.:1.0000  
##  Max.   :15.000           Max.   :13.000   Max.   :2.0000   Max.   :2.0000  
##                                                                             
##  total_embryos      prop_live     
##  Min.   : 6.000   Min.   :0.7500  
##  1st Qu.: 7.000   1st Qu.:0.8333  
##  Median : 9.000   Median :1.0000  
##  Mean   : 9.049   Mean   :0.9238  
##  3rd Qu.:10.000   3rd Qu.:1.0000  
##  Max.   :16.000   Max.   :1.0000  
## 

PS

All PS:

ps_all <- neonate_postshed_clean %>% 
  left_join(postshed_osmol,
            by = c('Indiv_ID', 'CEWL_Blood_Collect_Date')) %>% 
  left_join(all_data_CEWL, 
            by = c('Indiv_ID', 'CEWL_Blood_Collect_Date'))

ps_all_mother <- ps_all %>% 
  left_join(mother_values, by = "Mother_ID")

B vs P

subset of neonates/clutches with measurements at BOTH time points:

dab_ps_compare <- dab_all %>% 
  # only keep clutches with post-shed msmts
  filter(Mother_ID %in% ps_all$Mother_ID) %>%
  # remove some cols that do not have matches in ps data
  select(-prop_live, -total_embryos, -slugs, -dead_at_birth, -live_neonates, -Mother_Days_in_Treatment, -Mother_osml, -Mother_SVL, -Mother_mass, -Mother_CEWL) %>% 
  # add those post-shed data
  rbind(ps_all) %>% 
  # formatting
  mutate(timept = factor(timept))

summary(dab_ps_compare)
##     Indiv_ID  Date_Born_Shed       Tail_Length_cm      SVL_cm    
##  303    : 2   Min.   :2022-08-26   Min.   :1.200   Min.   :20.5  
##  304    : 2   1st Qu.:2022-08-28   1st Qu.:1.625   1st Qu.:22.3  
##  320    : 2   Median :2022-08-30   Median :1.950   Median :23.0  
##  300    : 1   Mean   :2022-08-31   Mean   :1.858   Mean   :23.1  
##  301    : 1   3rd Qu.:2022-09-03   3rd Qu.:2.000   3rd Qu.:24.0  
##  302    : 1   Max.   :2022-09-04   Max.   :2.500   Max.   :26.5  
##  (Other):65                                                      
##  Tail_SVL_Ratio        Mass_g      Sex      Mother_ID  Treatment
##  Min.   :0.05330   Min.   :11.10   F:24   125    :10   C:39     
##  1st Qu.:0.06945   1st Qu.:12.82   M:50   128    :10   W:35     
##  Median :0.08405   Median :14.15          133    :10            
##  Mean   :0.08059   Mean   :14.37          134    :10            
##  3rd Qu.:0.08930   3rd Qu.:15.28          135    :10            
##  Max.   :0.10330   Max.   :20.00          129    : 9            
##                                           (Other):15            
##    Tb_CEWL_c     Pee    CEWL_Blood_Collect_Date timept   Plasma_Osmol_Rep_mean
##  Min.   :25.80   N:54   Min.   :2022-08-28      DAB:39   Min.   :242.7        
##  1st Qu.:27.62   Y:20   1st Qu.:2022-08-29      PS :35   1st Qu.:286.6        
##  Median :28.50          Median :2022-08-31               Median :296.8        
##  Mean   :28.61          Mean   :2022-09-01               Mean   :304.3        
##  3rd Qu.:29.60          3rd Qu.:2022-09-04               3rd Qu.:315.9        
##  Max.   :31.30          Max.   :2022-09-05               Max.   :428.0        
##                                                                               
##    CEWL_g_m2h      msmt_temp_C    msmt_RH_percent      reps      
##  Min.   : 4.745   Min.   :23.56   Min.   :35.15   Min.   :3.000  
##  1st Qu.:10.897   1st Qu.:25.32   1st Qu.:38.77   1st Qu.:4.000  
##  Median :13.725   Median :25.73   Median :40.25   Median :4.000  
##  Mean   :13.829   Mean   :25.64   Mean   :40.51   Mean   :4.192  
##  3rd Qu.:16.712   3rd Qu.:26.10   3rd Qu.:42.20   3rd Qu.:4.000  
##  Max.   :23.808   Max.   :27.20   Max.   :46.22   Max.   :8.000  
##  NA's   :1        NA's   :1       NA's   :1       NA's   :1      
##     e_s_kPa         e_a_kPa         VPD_kPa     
##  Min.   :2.904   Min.   :1.148   Min.   :1.566  
##  1st Qu.:3.228   1st Qu.:1.258   1st Qu.:1.906  
##  Median :3.305   Median :1.330   Median :1.960  
##  Mean   :3.292   Mean   :1.332   Mean   :1.960  
##  3rd Qu.:3.380   3rd Qu.:1.400   3rd Qu.:2.059  
##  Max.   :3.606   Max.   :1.532   Max.   :2.179  
##  NA's   :1       NA's   :1       NA's   :1

Clutch Avgs

Get min/max, sd, and mean per clutch per variable.

DAB Proportions

dab_all_props <- dab_all %>% 
  #group_by(Mother_ID, Treatment, timept, Sex) %>%
  #count() %>% 
  group_by(Mother_ID, Treatment, timept) %>%
  summarise(
    n_neonates = n(), # number MEASURED, not n birthed
    Prop_Pee = sum(Pee == "Y")/n_neonates,
    Prop_F = sum(Sex == "F")/n_neonates
  )

summary(dab_all_props)
##    Mother_ID  Treatment    timept            n_neonates     Prop_Pee     
##  101    : 1   C: 8      Length:18          Min.   :3.0   Min.   :0.0000  
##  103    : 1   W:10      Class :character   1st Qu.:4.0   1st Qu.:0.2875  
##  104    : 1             Mode  :character   Median :5.0   Median :0.6000  
##  105    : 1                                Mean   :4.5   Mean   :0.5287  
##  112    : 1                                3rd Qu.:5.0   3rd Qu.:0.7875  
##  114    : 1                                Max.   :5.0   Max.   :1.0000  
##  (Other):12                                                              
##      Prop_F      
##  Min.   :0.0000  
##  1st Qu.:0.2000  
##  Median :0.3667  
##  Mean   :0.3296  
##  3rd Qu.:0.5000  
##  Max.   :0.7500  
## 

DAB Means

dab_all_clutch_avg <- dab_all %>% 
  group_by(Mother_ID, Treatment, timept) %>%
  mutate(n = n()) %>% 
  group_by(Mother_ID, Treatment, timept, n) %>% 
  summarise_at(
    vars(Tail_Length_cm, SVL_cm, Mass_g, Tb_CEWL_c, 
         Plasma_Osmol_Rep_mean, CEWL_g_m2h, 
         msmt_temp_C, msmt_RH_percent, VPD_kPa), 
    list(mean = mean, sd = sd, min = min, max = max),
    na.rm = T
  )

length(unique(dab_all_clutch_avg$Mother_ID)) == nrow(dab_all_clutch_avg)
## [1] TRUE
summary(dab_all_clutch_avg)
##    Mother_ID  Treatment    timept                n       Tail_Length_cm_mean
##  101    : 1   C: 8      Length:18          Min.   :3.0   Min.   :1.560      
##  103    : 1   W:10      Class :character   1st Qu.:4.0   1st Qu.:1.706      
##  104    : 1             Mode  :character   Median :5.0   Median :1.877      
##  105    : 1                                Mean   :4.5   Mean   :1.833      
##  112    : 1                                3rd Qu.:5.0   3rd Qu.:1.945      
##  114    : 1                                Max.   :5.0   Max.   :2.100      
##  (Other):12                                                                 
##   SVL_cm_mean     Mass_g_mean    Tb_CEWL_c_mean  Plasma_Osmol_Rep_mean_mean
##  Min.   :21.80   Min.   :12.80   Min.   :26.32   Min.   :280.5             
##  1st Qu.:22.08   1st Qu.:13.91   1st Qu.:27.76   1st Qu.:284.8             
##  Median :22.43   Median :14.92   Median :28.09   Median :297.8             
##  Mean   :22.67   Mean   :15.27   Mean   :28.24   Mean   :296.6             
##  3rd Qu.:23.07   3rd Qu.:16.06   3rd Qu.:28.73   3rd Qu.:302.7             
##  Max.   :24.94   Max.   :22.30   Max.   :29.92   Max.   :324.7             
##                                                                            
##  CEWL_g_m2h_mean  msmt_temp_C_mean msmt_RH_percent_mean  VPD_kPa_mean  
##  Min.   : 5.005   Min.   :23.71    Min.   :34.38        Min.   :1.595  
##  1st Qu.:10.705   1st Qu.:25.42    1st Qu.:38.78        1st Qu.:1.902  
##  Median :12.709   Median :26.03    Median :40.85        Median :1.959  
##  Mean   :12.076   Mean   :25.83    Mean   :40.68        Mean   :1.976  
##  3rd Qu.:14.229   3rd Qu.:26.33    3rd Qu.:41.42        3rd Qu.:2.094  
##  Max.   :15.892   Max.   :26.88    Max.   :46.73        Max.   :2.293  
##                                                                        
##  Tail_Length_cm_sd   SVL_cm_sd        Mass_g_sd       Tb_CEWL_c_sd   
##  Min.   :0.1633    Min.   :0.2708   Min.   :0.2380   Min.   :0.3768  
##  1st Qu.:0.2527    1st Qu.:0.4820   1st Qu.:0.3880   1st Qu.:0.5765  
##  Median :0.2791    Median :0.6630   Median :0.6510   Median :0.8608  
##  Mean   :0.2780    Mean   :0.6670   Mean   :0.8053   Mean   :0.9451  
##  3rd Qu.:0.3125    3rd Qu.:0.8281   3rd Qu.:1.0364   3rd Qu.:1.2762  
##  Max.   :0.4359    Max.   :1.0954   Max.   :2.0768   Max.   :1.8050  
##                                                                      
##  Plasma_Osmol_Rep_mean_sd CEWL_g_m2h_sd    msmt_temp_C_sd   msmt_RH_percent_sd
##  Min.   : 1.250           Min.   :0.7602   Min.   :0.1313   Min.   :0.05774   
##  1st Qu.: 6.050           1st Qu.:1.6437   1st Qu.:0.1575   1st Qu.:0.38341   
##  Median : 8.307           Median :2.0524   Median :0.1925   Median :0.58473   
##  Mean   :10.125           Mean   :2.4926   Mean   :0.2216   Mean   :0.82751   
##  3rd Qu.:11.687           3rd Qu.:2.8243   3rd Qu.:0.2910   3rd Qu.:1.03296   
##  Max.   :34.588           Max.   :5.4595   Max.   :0.3538   Max.   :2.60736   
##                                                                               
##    VPD_kPa_sd       Tail_Length_cm_min   SVL_cm_min      Mass_g_min   
##  Min.   :0.009671   Min.   :1.200      Min.   :20.50   Min.   :11.60  
##  1st Qu.:0.026563   1st Qu.:1.400      1st Qu.:21.32   1st Qu.:12.80  
##  Median :0.040532   Median :1.500      Median :21.60   Median :13.40  
##  Mean   :0.043391   Mean   :1.506      Mean   :21.79   Mean   :14.26  
##  3rd Qu.:0.056447   3rd Qu.:1.600      3rd Qu.:22.43   3rd Qu.:14.90  
##  Max.   :0.110722   Max.   :1.900      Max.   :24.00   Max.   :21.50  
##                                                                       
##  Tb_CEWL_c_min   Plasma_Osmol_Rep_mean_min CEWL_g_m2h_min   msmt_temp_C_min
##  Min.   :25.80   Min.   :263.0             Min.   : 4.043   Min.   :23.56  
##  1st Qu.:26.82   1st Qu.:280.4             1st Qu.: 8.108   1st Qu.:25.24  
##  Median :27.05   Median :285.8             Median : 9.804   Median :25.82  
##  Mean   :27.22   Mean   :284.8             Mean   : 9.534   Mean   :25.57  
##  3rd Qu.:27.65   3rd Qu.:288.9             3rd Qu.:10.988   3rd Qu.:26.06  
##  Max.   :29.50   Max.   :302.5             Max.   :12.994   Max.   :26.70  
##                                                                            
##  msmt_RH_percent_min  VPD_kPa_min    Tail_Length_cm_max   SVL_cm_max   
##  Min.   :32.92       Min.   :1.566   Min.   :2.000      Min.   :22.40  
##  1st Qu.:37.81       1st Qu.:1.837   1st Qu.:2.000      1st Qu.:23.00  
##  Median :40.20       Median :1.918   Median :2.050      Median :23.30  
##  Mean   :39.78       Mean   :1.923   Mean   :2.122      Mean   :23.46  
##  3rd Qu.:41.22       3rd Qu.:2.040   3rd Qu.:2.275      3rd Qu.:23.80  
##  Max.   :45.88       Max.   :2.214   Max.   :2.500      Max.   :25.70  
##                                                                        
##    Mass_g_max    Tb_CEWL_c_max   Plasma_Osmol_Rep_mean_max CEWL_g_m2h_max  
##  Min.   :13.30   Min.   :27.30   Min.   :282.0             Min.   : 5.973  
##  1st Qu.:14.53   1st Qu.:28.62   1st Qu.:293.1             1st Qu.:12.828  
##  Median :16.00   Median :29.45   Median :307.2             Median :15.102  
##  Mean   :16.13   Mean   :29.46   Mean   :307.9             Mean   :15.398  
##  3rd Qu.:17.05   3rd Qu.:30.27   3rd Qu.:315.5             3rd Qu.:17.650  
##  Max.   :23.00   Max.   :31.30   Max.   :355.0             Max.   :23.808  
##                                                                            
##  msmt_temp_C_max msmt_RH_percent_max  VPD_kPa_max   
##  Min.   :24.02   Min.   :36.30       Min.   :1.641  
##  1st Qu.:25.64   1st Qu.:39.56       1st Qu.:1.933  
##  Median :26.29   Median :41.69       Median :2.004  
##  Mean   :26.09   Mean   :41.71       Mean   :2.024  
##  3rd Qu.:26.59   3rd Qu.:42.69       3rd Qu.:2.143  
##  Max.   :27.20   Max.   :49.42       Max.   :2.324  
## 

Add Mother Data!

dab_clutch_mom <- dab_all_clutch_avg %>% 
  left_join(mother_values, by = "Mother_ID") %>% 
  left_join(dab_all_props, by = c("Mother_ID", "Treatment", "timept")) %>% 
  left_join(neo_counts, by = "Mother_ID") %>% 
  mutate(Mother_mass_kg = Mother_mass/1000)

summary(dab_clutch_mom)
##    Mother_ID  Treatment    timept                n       Tail_Length_cm_mean
##  101    : 1   C: 8      Length:18          Min.   :3.0   Min.   :1.560      
##  103    : 1   W:10      Class :character   1st Qu.:4.0   1st Qu.:1.706      
##  104    : 1             Mode  :character   Median :5.0   Median :1.877      
##  105    : 1                                Mean   :4.5   Mean   :1.833      
##  112    : 1                                3rd Qu.:5.0   3rd Qu.:1.945      
##  114    : 1                                Max.   :5.0   Max.   :2.100      
##  (Other):12                                                                 
##   SVL_cm_mean     Mass_g_mean    Tb_CEWL_c_mean  Plasma_Osmol_Rep_mean_mean
##  Min.   :21.80   Min.   :12.80   Min.   :26.32   Min.   :280.5             
##  1st Qu.:22.08   1st Qu.:13.91   1st Qu.:27.76   1st Qu.:284.8             
##  Median :22.43   Median :14.92   Median :28.09   Median :297.8             
##  Mean   :22.67   Mean   :15.27   Mean   :28.24   Mean   :296.6             
##  3rd Qu.:23.07   3rd Qu.:16.06   3rd Qu.:28.73   3rd Qu.:302.7             
##  Max.   :24.94   Max.   :22.30   Max.   :29.92   Max.   :324.7             
##                                                                            
##  CEWL_g_m2h_mean  msmt_temp_C_mean msmt_RH_percent_mean  VPD_kPa_mean  
##  Min.   : 5.005   Min.   :23.71    Min.   :34.38        Min.   :1.595  
##  1st Qu.:10.705   1st Qu.:25.42    1st Qu.:38.78        1st Qu.:1.902  
##  Median :12.709   Median :26.03    Median :40.85        Median :1.959  
##  Mean   :12.076   Mean   :25.83    Mean   :40.68        Mean   :1.976  
##  3rd Qu.:14.229   3rd Qu.:26.33    3rd Qu.:41.42        3rd Qu.:2.094  
##  Max.   :15.892   Max.   :26.88    Max.   :46.73        Max.   :2.293  
##                                                                        
##  Tail_Length_cm_sd   SVL_cm_sd        Mass_g_sd       Tb_CEWL_c_sd   
##  Min.   :0.1633    Min.   :0.2708   Min.   :0.2380   Min.   :0.3768  
##  1st Qu.:0.2527    1st Qu.:0.4820   1st Qu.:0.3880   1st Qu.:0.5765  
##  Median :0.2791    Median :0.6630   Median :0.6510   Median :0.8608  
##  Mean   :0.2780    Mean   :0.6670   Mean   :0.8053   Mean   :0.9451  
##  3rd Qu.:0.3125    3rd Qu.:0.8281   3rd Qu.:1.0364   3rd Qu.:1.2762  
##  Max.   :0.4359    Max.   :1.0954   Max.   :2.0768   Max.   :1.8050  
##                                                                      
##  Plasma_Osmol_Rep_mean_sd CEWL_g_m2h_sd    msmt_temp_C_sd   msmt_RH_percent_sd
##  Min.   : 1.250           Min.   :0.7602   Min.   :0.1313   Min.   :0.05774   
##  1st Qu.: 6.050           1st Qu.:1.6437   1st Qu.:0.1575   1st Qu.:0.38341   
##  Median : 8.307           Median :2.0524   Median :0.1925   Median :0.58473   
##  Mean   :10.125           Mean   :2.4926   Mean   :0.2216   Mean   :0.82751   
##  3rd Qu.:11.687           3rd Qu.:2.8243   3rd Qu.:0.2910   3rd Qu.:1.03296   
##  Max.   :34.588           Max.   :5.4595   Max.   :0.3538   Max.   :2.60736   
##                                                                               
##    VPD_kPa_sd       Tail_Length_cm_min   SVL_cm_min      Mass_g_min   
##  Min.   :0.009671   Min.   :1.200      Min.   :20.50   Min.   :11.60  
##  1st Qu.:0.026563   1st Qu.:1.400      1st Qu.:21.32   1st Qu.:12.80  
##  Median :0.040532   Median :1.500      Median :21.60   Median :13.40  
##  Mean   :0.043391   Mean   :1.506      Mean   :21.79   Mean   :14.26  
##  3rd Qu.:0.056447   3rd Qu.:1.600      3rd Qu.:22.43   3rd Qu.:14.90  
##  Max.   :0.110722   Max.   :1.900      Max.   :24.00   Max.   :21.50  
##                                                                       
##  Tb_CEWL_c_min   Plasma_Osmol_Rep_mean_min CEWL_g_m2h_min   msmt_temp_C_min
##  Min.   :25.80   Min.   :263.0             Min.   : 4.043   Min.   :23.56  
##  1st Qu.:26.82   1st Qu.:280.4             1st Qu.: 8.108   1st Qu.:25.24  
##  Median :27.05   Median :285.8             Median : 9.804   Median :25.82  
##  Mean   :27.22   Mean   :284.8             Mean   : 9.534   Mean   :25.57  
##  3rd Qu.:27.65   3rd Qu.:288.9             3rd Qu.:10.988   3rd Qu.:26.06  
##  Max.   :29.50   Max.   :302.5             Max.   :12.994   Max.   :26.70  
##                                                                            
##  msmt_RH_percent_min  VPD_kPa_min    Tail_Length_cm_max   SVL_cm_max   
##  Min.   :32.92       Min.   :1.566   Min.   :2.000      Min.   :22.40  
##  1st Qu.:37.81       1st Qu.:1.837   1st Qu.:2.000      1st Qu.:23.00  
##  Median :40.20       Median :1.918   Median :2.050      Median :23.30  
##  Mean   :39.78       Mean   :1.923   Mean   :2.122      Mean   :23.46  
##  3rd Qu.:41.22       3rd Qu.:2.040   3rd Qu.:2.275      3rd Qu.:23.80  
##  Max.   :45.88       Max.   :2.214   Max.   :2.500      Max.   :25.70  
##                                                                        
##    Mass_g_max    Tb_CEWL_c_max   Plasma_Osmol_Rep_mean_max CEWL_g_m2h_max  
##  Min.   :13.30   Min.   :27.30   Min.   :282.0             Min.   : 5.973  
##  1st Qu.:14.53   1st Qu.:28.62   1st Qu.:293.1             1st Qu.:12.828  
##  Median :16.00   Median :29.45   Median :307.2             Median :15.102  
##  Mean   :16.13   Mean   :29.46   Mean   :307.9             Mean   :15.398  
##  3rd Qu.:17.05   3rd Qu.:30.27   3rd Qu.:315.5             3rd Qu.:17.650  
##  Max.   :23.00   Max.   :31.30   Max.   :355.0             Max.   :23.808  
##                                                                            
##  msmt_temp_C_max msmt_RH_percent_max  VPD_kPa_max     Mother_CEWL   
##  Min.   :24.02   Min.   :36.30       Min.   :1.641   Min.   :10.27  
##  1st Qu.:25.64   1st Qu.:39.56       1st Qu.:1.933   1st Qu.:17.18  
##  Median :26.29   Median :41.69       Median :2.004   Median :19.82  
##  Mean   :26.09   Mean   :41.71       Mean   :2.024   Mean   :19.22  
##  3rd Qu.:26.59   3rd Qu.:42.69       3rd Qu.:2.143   3rd Qu.:21.22  
##  Max.   :27.20   Max.   :49.42       Max.   :2.324   Max.   :27.82  
##                                                                     
##   Mother_mass      Mother_SVL     Mother_osml    Mother_Days_in_Treatment
##  Min.   :181.7   Min.   :72.00   Min.   :281.7   Min.   : 2.000          
##  1st Qu.:281.4   1st Qu.:75.75   1st Qu.:304.2   1st Qu.: 5.000          
##  Median :314.3   Median :80.90   Median :312.0   Median : 8.000          
##  Mean   :334.6   Mean   :81.32   Mean   :313.5   Mean   : 8.833          
##  3rd Qu.:371.0   3rd Qu.:85.88   3rd Qu.:320.6   3rd Qu.:14.000          
##  Max.   :570.9   Max.   :91.00   Max.   :357.2   Max.   :15.000          
##                                                                          
##    n_neonates     Prop_Pee          Prop_F       live_neonates   
##  Min.   :3.0   Min.   :0.0000   Min.   :0.0000   Min.   : 5.000  
##  1st Qu.:4.0   1st Qu.:0.2875   1st Qu.:0.2000   1st Qu.: 7.000  
##  Median :5.0   Median :0.6000   Median :0.3667   Median : 8.000  
##  Mean   :4.5   Mean   :0.5287   Mean   :0.3296   Mean   : 8.278  
##  3rd Qu.:5.0   3rd Qu.:0.7875   3rd Qu.:0.5000   3rd Qu.: 9.750  
##  Max.   :5.0   Max.   :1.0000   Max.   :0.7500   Max.   :13.000  
##                                                                  
##  dead_at_birth        slugs        total_embryos     prop_live     
##  Min.   :0.0000   Min.   :0.0000   Min.   : 6.00   Min.   :0.7500  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.: 7.25   1st Qu.:0.8472  
##  Median :0.0000   Median :0.0000   Median : 9.00   Median :1.0000  
##  Mean   :0.3889   Mean   :0.3333   Mean   : 9.00   Mean   :0.9268  
##  3rd Qu.:0.7500   3rd Qu.:0.7500   3rd Qu.:10.00   3rd Qu.:1.0000  
##  Max.   :2.0000   Max.   :2.0000   Max.   :16.00   Max.   :1.0000  
##                                                                    
##  Mother_mass_kg  
##  Min.   :0.1817  
##  1st Qu.:0.2814  
##  Median :0.3143  
##  Mean   :0.3346  
##  3rd Qu.:0.3710  
##  Max.   :0.5709  
## 

PS Means

ps_all_clutch_avg <- ps_all %>% 
  group_by(Mother_ID, Treatment, timept) %>%
  mutate(n = n()) %>% 
  group_by(Mother_ID, Treatment, timept, n) %>% 
  summarise_at(
    vars(Tail_Length_cm, SVL_cm, Mass_g, Tb_CEWL_c, 
         Plasma_Osmol_Rep_mean, CEWL_g_m2h, 
         msmt_temp_C, msmt_RH_percent, VPD_kPa), 
    list(mean = mean, sd = sd, min = min, max = max),
    na.rm = T
  )

length(unique(ps_all_clutch_avg$Mother_ID)) == nrow(ps_all_clutch_avg)
## [1] TRUE

B vs P

dab_ps_compare_means <- dab_all_clutch_avg %>% 
  # only keep clutches with post-shed msmts
  filter(Mother_ID %in% ps_all$Mother_ID) %>%
  # add those post-shed data
  rbind(ps_all_clutch_avg) %>% 
  # formatting
  mutate(timept = factor(timept))

summary(dab_ps_compare_means)
##    Mother_ID Treatment timept        n         Tail_Length_cm_mean
##  124    :2   C:8       DAB:8   Min.   :3.000   Min.   :1.560      
##  125    :2   W:8       PS :8   1st Qu.:4.750   1st Qu.:1.775      
##  128    :2                     Median :5.000   Median :1.880      
##  129    :2                     Mean   :4.625   Mean   :1.865      
##  131    :2                     3rd Qu.:5.000   3rd Qu.:1.962      
##  133    :2                     Max.   :5.000   Max.   :2.160      
##  (Other):4                                                        
##   SVL_cm_mean     Mass_g_mean    Tb_CEWL_c_mean  Plasma_Osmol_Rep_mean_mean
##  Min.   :21.80   Min.   :11.66   Min.   :27.42   Min.   :280.5             
##  1st Qu.:22.31   1st Qu.:13.07   1st Qu.:27.96   1st Qu.:286.2             
##  Median :23.35   Median :14.14   Median :28.60   Median :296.6             
##  Mean   :23.17   Mean   :14.40   Mean   :28.58   Mean   :304.2             
##  3rd Qu.:23.72   3rd Qu.:15.07   3rd Qu.:29.23   3rd Qu.:312.5             
##  Max.   :25.87   Max.   :18.46   Max.   :29.92   Max.   :363.4             
##                                                                            
##  CEWL_g_m2h_mean  msmt_temp_C_mean msmt_RH_percent_mean  VPD_kPa_mean  
##  Min.   : 7.296   Min.   :23.71    Min.   :36.80        Min.   :1.595  
##  1st Qu.:10.795   1st Qu.:25.40    1st Qu.:39.02        1st Qu.:1.923  
##  Median :15.046   Median :25.75    Median :40.10        Median :1.962  
##  Mean   :13.959   Mean   :25.65    Mean   :40.39        Mean   :1.965  
##  3rd Qu.:16.193   3rd Qu.:26.01    3rd Qu.:41.24        3rd Qu.:2.060  
##  Max.   :17.244   Max.   :26.88    Max.   :45.57        Max.   :2.160  
##                                                                        
##  Tail_Length_cm_sd   SVL_cm_sd        Mass_g_sd       Tb_CEWL_c_sd   
##  Min.   :0.1517    Min.   :0.4435   Min.   :0.2380   Min.   :0.3000  
##  1st Qu.:0.1930    1st Qu.:0.6132   1st Qu.:0.4584   1st Qu.:0.4959  
##  Median :0.2334    Median :0.7606   Median :0.6737   Median :0.9314  
##  Mean   :0.2407    Mean   :0.7866   Mean   :0.8755   Mean   :0.9654  
##  3rd Qu.:0.2842    3rd Qu.:0.9680   3rd Qu.:1.2115   3rd Qu.:1.3562  
##  Max.   :0.3808    Max.   :1.0970   Max.   :2.0768   Max.   :1.8050  
##                                                                      
##  Plasma_Osmol_Rep_mean_sd CEWL_g_m2h_sd   msmt_temp_C_sd   msmt_RH_percent_sd
##  Min.   : 2.739           Min.   :1.128   Min.   :0.1311   Min.   :0.1248    
##  1st Qu.: 6.355           1st Qu.:1.765   1st Qu.:0.1642   1st Qu.:0.3391    
##  Median : 8.009           Median :2.213   Median :0.1990   Median :0.6041    
##  Mean   :12.556           Mean   :2.577   Mean   :0.2269   Mean   :0.8738    
##  3rd Qu.:14.390           3rd Qu.:3.172   3rd Qu.:0.2876   3rd Qu.:1.2300    
##  Max.   :43.634           Max.   :5.040   Max.   :0.3929   Max.   :2.6074    
##                                                                              
##    VPD_kPa_sd       Tail_Length_cm_min   SVL_cm_min      Mass_g_min   
##  Min.   :0.009671   Min.   :1.200      Min.   :20.50   Min.   :11.10  
##  1st Qu.:0.028434   1st Qu.:1.500      1st Qu.:21.38   1st Qu.:11.90  
##  Median :0.037982   Median :1.500      Median :22.45   Median :13.05  
##  Mean   :0.044378   Mean   :1.562      Mean   :22.18   Mean   :13.23  
##  3rd Qu.:0.062645   3rd Qu.:1.700      3rd Qu.:23.00   3rd Qu.:14.28  
##  Max.   :0.110722   Max.   :1.900      Max.   :24.60   Max.   :16.60  
##                                                                       
##  Tb_CEWL_c_min   Plasma_Osmol_Rep_mean_min CEWL_g_m2h_min   msmt_temp_C_min
##  Min.   :25.80   Min.   :242.7             Min.   : 4.745   Min.   :23.56  
##  1st Qu.:26.98   1st Qu.:279.5             1st Qu.: 9.571   1st Qu.:25.06  
##  Median :27.40   Median :286.2             Median :11.494   Median :25.52  
##  Mean   :27.50   Mean   :290.2             Mean   :10.959   Mean   :25.40  
##  3rd Qu.:27.93   3rd Qu.:296.2             3rd Qu.:12.982   3rd Qu.:25.76  
##  Max.   :29.50   Max.   :332.5             Max.   :14.812   Max.   :26.70  
##                                                                            
##  msmt_RH_percent_min  VPD_kPa_min    Tail_Length_cm_max   SVL_cm_max   
##  Min.   :35.15       Min.   :1.566   Min.   :2.000      Min.   :22.40  
##  1st Qu.:37.66       1st Qu.:1.858   1st Qu.:2.000      1st Qu.:23.45  
##  Median :39.23       Median :1.917   Median :2.050      Median :24.00  
##  Mean   :39.44       Mean   :1.911   Mean   :2.119      Mean   :24.11  
##  3rd Qu.:40.96       3rd Qu.:1.972   3rd Qu.:2.200      3rd Qu.:24.62  
##  Max.   :44.34       Max.   :2.144   Max.   :2.500      Max.   :26.50  
##                                                                        
##    Mass_g_max    Tb_CEWL_c_max   Plasma_Osmol_Rep_mean_max CEWL_g_m2h_max  
##  Min.   :12.10   Min.   :28.10   Min.   :283.5             Min.   : 9.705  
##  1st Qu.:13.70   1st Qu.:29.35   1st Qu.:297.5             1st Qu.:13.229  
##  Median :14.70   Median :30.20   Median :306.2             Median :18.466  
##  Mean   :15.25   Mean   :29.83   Mean   :318.5             Mean   :17.158  
##  3rd Qu.:16.38   3rd Qu.:30.43   3rd Qu.:327.8             3rd Qu.:19.771  
##  Max.   :20.00   Max.   :31.30   Max.   :428.0             Max.   :23.808  
##                                                                            
##  msmt_temp_C_max msmt_RH_percent_max  VPD_kPa_max   
##  Min.   :24.02   Min.   :37.77       Min.   :1.641  
##  1st Qu.:25.68   1st Qu.:40.20       1st Qu.:1.963  
##  Median :25.99   Median :41.31       Median :2.001  
##  Mean   :25.93   Mean   :41.46       Mean   :2.014  
##  3rd Qu.:26.29   3rd Qu.:42.48       3rd Qu.:2.138  
##  Max.   :27.20   Max.   :46.22       Max.   :2.179  
## 

Exp Info

dab_ps_compare %>% 
  group_by(Mother_ID) %>% 
  summarise(
    min_date = min(CEWL_Blood_Collect_Date),
    max_date = max(CEWL_Blood_Collect_Date)
  ) %>% 
  mutate(range = max_date - min_date) %>% 
  arrange(range)
## # A tibble: 8 × 4
##   Mother_ID min_date   max_date   range 
##   <fct>     <date>     <date>     <drtn>
## 1 124       2022-08-31 2022-09-05 5 days
## 2 131       2022-08-31 2022-09-05 5 days
## 3 133       2022-08-31 2022-09-05 5 days
## 4 134       2022-08-31 2022-09-05 5 days
## 5 128       2022-08-29 2022-09-04 6 days
## 6 129       2022-08-28 2022-09-03 6 days
## 7 135       2022-08-28 2022-09-03 6 days
## 8 125       2022-08-28 2022-09-04 7 days
dab_ps_compare %>% 
  summarise(
    min_date = min(CEWL_Blood_Collect_Date),
    max_date = max(CEWL_Blood_Collect_Date)
  ) %>% 
  mutate(range = max_date - min_date)
##     min_date   max_date  range
## 1 2022-08-28 2022-09-05 8 days

Test Repro Output

summary(lm(data = dab_clutch_mom, Prop_Pee~Treatment))
## 
## Call:
## lm(formula = Prop_Pee ~ Treatment, data = dab_clutch_mom)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4500 -0.3109  0.0750  0.2604  0.4354 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   0.3146     0.1027   3.063  0.00744 **
## TreatmentW    0.3854     0.1378   2.797  0.01292 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2905 on 16 degrees of freedom
## Multiple R-squared:  0.3284, Adjusted R-squared:  0.2864 
## F-statistic: 7.823 on 1 and 16 DF,  p-value: 0.01292
summary(lm(data = dab_clutch_mom, Prop_F~Treatment))
## 
## Call:
## lm(formula = Prop_F ~ Treatment, data = dab_clutch_mom)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.39167 -0.16375 -0.01083  0.18625  0.47000 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.39167    0.08315   4.711 0.000236 ***
## TreatmentW  -0.11167    0.11155  -1.001 0.331711    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2352 on 16 degrees of freedom
## Multiple R-squared:  0.05894,    Adjusted R-squared:  0.0001214 
## F-statistic: 1.002 on 1 and 16 DF,  p-value: 0.3317
summary(lm(data = dab_clutch_mom, live_neonates~Treatment))
## 
## Call:
## lm(formula = live_neonates ~ Treatment, data = dab_clutch_mom)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -3.10  -1.50  -0.10   1.25   4.50 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    8.500      0.765   11.11 6.22e-09 ***
## TreatmentW    -0.400      1.026   -0.39    0.702    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.164 on 16 degrees of freedom
## Multiple R-squared:  0.009405,   Adjusted R-squared:  -0.05251 
## F-statistic: 0.1519 on 1 and 16 DF,  p-value: 0.7019
summary(lm(data = dab_clutch_mom, dead_at_birth~Treatment))
## 
## Call:
## lm(formula = dead_at_birth ~ Treatment, data = dab_clutch_mom)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -0.50  -0.45  -0.30   0.30   1.70 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   0.5000     0.2516   1.988   0.0642 .
## TreatmentW   -0.2000     0.3375  -0.593   0.5617  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7115 on 16 degrees of freedom
## Multiple R-squared:  0.02148,    Adjusted R-squared:  -0.03968 
## F-statistic: 0.3512 on 1 and 16 DF,  p-value: 0.5617
summary(lm(data = dab_clutch_mom, slugs~Treatment))
## 
## Call:
## lm(formula = slugs ~ Treatment, data = dab_clutch_mom)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -0.500 -0.425 -0.200  0.325  1.500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   0.5000     0.2092   2.390   0.0295 *
## TreatmentW   -0.3000     0.2806  -1.069   0.3009  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5916 on 16 degrees of freedom
## Multiple R-squared:  0.06667,    Adjusted R-squared:  0.008333 
## F-statistic: 1.143 on 1 and 16 DF,  p-value: 0.3009
summary(lm(data = dab_clutch_mom, total_embryos~Treatment))
## 
## Call:
## lm(formula = total_embryos ~ Treatment, data = dab_clutch_mom)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.600 -2.250 -0.500  1.175  6.500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.5000     0.8856  10.727 1.03e-08 ***
## TreatmentW   -0.9000     1.1882  -0.757     0.46    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.505 on 16 degrees of freedom
## Multiple R-squared:  0.03462,    Adjusted R-squared:  -0.02572 
## F-statistic: 0.5737 on 1 and 16 DF,  p-value: 0.4598
still_neo_mom_mass_lm <- lm(data = dab_clutch_mom, 
                  dead_at_birth ~ Treatment * Mother_mass)
car::Anova(still_neo_mom_mass_lm, type = "2", test = "F")
## Anova Table (Type II tests)
## 
## Response: dead_at_birth
##                       Sum Sq Df F value Pr(>F)
## Treatment             0.3236  1  0.5845 0.4573
## Mother_mass           0.3477  1  0.6281 0.4413
## Treatment:Mother_mass 0.0016  1  0.0028 0.9585
## Residuals             7.7507 14
slug_mom_mass_lm <- lm(data = dab_clutch_mom, 
                  slugs ~ Treatment * Mother_mass)
car::Anova(slug_mom_mass_lm, type = "2", test = "F")
## Anova Table (Type II tests)
## 
## Response: slugs
##                       Sum Sq Df F value Pr(>F)
## Treatment             0.3003  1  0.8641 0.3684
## Mother_mass           0.0459  1  0.1322 0.7216
## Treatment:Mother_mass 0.6880  1  1.9794 0.1813
## Residuals             4.8660 14

Live ~ Tmt * Mass

live_neo_mom_mass_lm <- lm(data = dab_clutch_mom, 
                  live_neonates ~ Treatment * Mother_mass)

summary(live_neo_mom_mass_lm)
## 
## Call:
## lm(formula = live_neonates ~ Treatment * Mother_mass, data = dab_clutch_mom)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9782 -0.7241  0.0131  0.6161  3.9702 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)   
## (Intercept)             0.59933    2.63951   0.227  0.82366   
## TreatmentW              6.85106    3.32457   2.061  0.05842 . 
## Mother_mass             0.02592    0.00841   3.082  0.00811 **
## TreatmentW:Mother_mass -0.02411    0.01000  -2.410  0.03027 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.781 on 14 degrees of freedom
## Multiple R-squared:  0.4127, Adjusted R-squared:  0.2868 
## F-statistic: 3.279 on 3 and 14 DF,  p-value: 0.05271
car::Anova(live_neo_mom_mass_lm, type = "2", test = "F")
## Anova Table (Type II tests)
## 
## Response: live_neonates
##                       Sum Sq Df F value  Pr(>F)  
## Treatment              3.150  1  0.9931 0.33590  
## Mother_mass           12.064  1  3.8031 0.07147 .
## Treatment:Mother_mass 18.427  1  5.8090 0.03027 *
## Residuals             44.409 14                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emt <- emtrends(live_neo_mom_mass_lm, "Treatment", var = "Mother_mass")
tidy(emt)
## # A tibble: 2 × 6
##   Treatment Mother_mass.trend std.error    df statistic p.value
##   <chr>                 <dbl>     <dbl> <dbl>     <dbl>   <dbl>
## 1 C                   0.0259    0.00841    14     3.08  0.00811
## 2 W                   0.00181   0.00542    14     0.335 0.743
pairs(emt)
##  contrast estimate   SE df t.ratio p.value
##  C - W      0.0241 0.01 14   2.410  0.0303
0.025922949*100 # change in mother grams
## [1] 2.592295
0.001812418*100
## [1] 0.1812418

Plot

plot(effect("Treatment:Mother_mass", live_neo_mom_mass_lm), multiline = T)

ggplot(
  dab_clutch_mom,
  aes(x = Mother_mass, y = live_neonates, color = Treatment)
) +
  geom_point() +
  geom_text(
    aes(label = Mother_ID),
    nudge_x = 0.25, nudge_y = 0.25, 
    check_overlap = T, show.legend = F
  ) +
  geom_smooth(
    method = "lm",
    se = F,
    formula = 'y~x'
  )

neo_counts %>% 
  filter(slugs > 0 | dead_at_birth > 0)
## # A tibble: 9 × 6
##   Mother_ID live_neonates dead_at_birth slugs total_embryos prop_live
##   <fct>             <dbl>         <dbl> <dbl>         <dbl>     <dbl>
## 1 133                  13             1     2            16     0.812
## 2 126                  10             1     0            11     0.909
## 3 125                   8             2     0            10     0.8  
## 4 128                   8             2     0            10     0.8  
## 5 129                   8             0     1             9     0.889
## 6 134                   8             0     1             9     0.889
## 7 105                   6             1     1             8     0.75 
## 8 117                   5             2     1             8     0.625
## 9 131                   5             0     1             6     0.833
prop_live_mom_mass_lm <- lm(data = dab_clutch_mom, 
                  prop_live ~ Treatment * Mother_mass)
summary(prop_live_mom_mass_lm)
## 
## Call:
## lm(formula = prop_live ~ Treatment * Mother_mass, data = dab_clutch_mom)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14775 -0.09115  0.04526  0.06506  0.09820 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            8.850e-01  1.453e-01   6.092 2.78e-05 ***
## TreatmentW             1.954e-02  1.830e-01   0.107    0.916    
## Mother_mass            7.043e-05  4.629e-04   0.152    0.881    
## TreatmentW:Mother_mass 3.736e-05  5.506e-04   0.068    0.947    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09802 on 14 degrees of freedom
## Multiple R-squared:  0.05259,    Adjusted R-squared:  -0.1504 
## F-statistic: 0.2591 on 3 and 14 DF,  p-value: 0.8537
car::Anova(prop_live_mom_mass_lm, type = "2", test = "F")
## Anova Table (Type II tests)
## 
## Response: prop_live
##                         Sum Sq Df F value Pr(>F)
## Treatment             0.004073  1  0.4239 0.5255
## Mother_mass           0.001435  1  0.1493 0.7050
## Treatment:Mother_mass 0.000044  1  0.0046 0.9469
## Residuals             0.134518 14

Live ~ SVL

live_neo_mom_SVL_lm <- lm(data = dab_clutch_mom, 
                  live_neonates ~ Mother_SVL)
summary(live_neo_mom_SVL_lm)
## 
## Call:
## lm(formula = live_neonates ~ Mother_SVL, data = dab_clutch_mom)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1395 -1.1640  0.0267  1.0432  3.4091 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -5.24127    5.95708  -0.880    0.392  
## Mother_SVL   0.16625    0.07305   2.276    0.037 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.889 on 16 degrees of freedom
## Multiple R-squared:  0.2445, Adjusted R-squared:  0.1973 
## F-statistic: 5.179 on 1 and 16 DF,  p-value: 0.03696
car::Anova(live_neo_mom_SVL_lm, type = "2", test = "F")
## Anova Table (Type II tests)
## 
## Response: live_neonates
##            Sum Sq Df F value  Pr(>F)  
## Mother_SVL 18.490  1  5.1792 0.03696 *
## Residuals  57.121 16                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
prop_live_mom_SVL_lm <- lm(data = dab_clutch_mom, 
                  prop_live ~ Mother_SVL)
summary(prop_live_mom_SVL_lm)
## 
## Call:
## lm(formula = prop_live ~ Mother_SVL, data = dab_clutch_mom)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.16830 -0.08855  0.06075  0.07283  0.08828 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  1.060478   0.295104   3.594  0.00243 **
## Mother_SVL  -0.001644   0.003619  -0.454  0.65579   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0936 on 16 degrees of freedom
## Multiple R-squared:  0.01273,    Adjusted R-squared:  -0.04897 
## F-statistic: 0.2063 on 1 and 16 DF,  p-value: 0.6558
car::Anova(prop_live_mom_SVL_lm, type = "2", test = "F")
## Anova Table (Type II tests)
## 
## Response: prop_live
##              Sum Sq Df F value Pr(>F)
## Mother_SVL 0.001807  1  0.2063 0.6558
## Residuals  0.140178 16
dab_clutch_mom %>% 
  group_by(Treatment, slugs, dead_at_birth) %>% 
  count()
## # A tibble: 9 × 4
## # Groups:   Treatment, slugs, dead_at_birth [9]
##   Treatment slugs dead_at_birth     n
##   <fct>     <dbl>         <dbl> <int>
## 1 C             0             0     4
## 2 C             0             2     1
## 3 C             1             0     1
## 4 C             1             1     1
## 5 C             2             1     1
## 6 W             0             0     6
## 7 W             0             1     1
## 8 W             0             2     1
## 9 W             1             0     2

DAB CEWL clutch ~ mother

this is what was done originally, but we now use all neonate data and do not avg by clutch

DAB_CEWL_LM_1 <- lm(data = dab_clutch_mom,
                      CEWL_g_m2h_mean ~ 
                      Treatment *
                      Mother_CEWL +
                      Tb_CEWL_c_mean +
                      msmt_temp_C_mean +
                      msmt_RH_percent_mean + 
                      VPD_kPa_mean)
                     
car::vif(DAB_CEWL_LM_1, type = "predictor")
##                             GVIF Df GVIF^(1/(2*Df)) Interacts With
## Treatment               4.744030  3        1.296257    Mother_CEWL
## Mother_CEWL             4.744030  3        1.296257      Treatment
## Tb_CEWL_c_mean          1.590281  1        1.261064           --  
## msmt_temp_C_mean      412.368805  1       20.306866           --  
## msmt_RH_percent_mean  454.046656  1       21.308371           --  
## VPD_kPa_mean         1121.024341  1       33.481702           --  
##                                                                                    Other Predictors
## Treatment                      Tb_CEWL_c_mean, msmt_temp_C_mean, msmt_RH_percent_mean, VPD_kPa_mean
## Mother_CEWL                    Tb_CEWL_c_mean, msmt_temp_C_mean, msmt_RH_percent_mean, VPD_kPa_mean
## Tb_CEWL_c_mean         Treatment, Mother_CEWL, msmt_temp_C_mean, msmt_RH_percent_mean, VPD_kPa_mean
## msmt_temp_C_mean         Treatment, Mother_CEWL, Tb_CEWL_c_mean, msmt_RH_percent_mean, VPD_kPa_mean
## msmt_RH_percent_mean         Treatment, Mother_CEWL, Tb_CEWL_c_mean, msmt_temp_C_mean, VPD_kPa_mean
## VPD_kPa_mean         Treatment, Mother_CEWL, Tb_CEWL_c_mean, msmt_temp_C_mean, msmt_RH_percent_mean
drop1(DAB_CEWL_LM_1)
## Single term deletions
## 
## Model:
## CEWL_g_m2h_mean ~ Treatment * Mother_CEWL + Tb_CEWL_c_mean + 
##     msmt_temp_C_mean + msmt_RH_percent_mean + VPD_kPa_mean
##                       Df Sum of Sq    RSS    AIC
## <none>                             46.632 33.134
## Tb_CEWL_c_mean         1    1.5292 48.161 31.715
## msmt_temp_C_mean       1    5.5176 52.149 33.147
## msmt_RH_percent_mean   1    5.8505 52.482 33.262
## VPD_kPa_mean           1    5.7394 52.371 33.224
## Treatment:Mother_CEWL  1    3.5723 50.204 32.463
car::Anova(DAB_CEWL_LM_1, type = "2")
## Anova Table (Type II tests)
## 
## Response: CEWL_g_m2h_mean
##                       Sum Sq Df F value Pr(>F)
## Treatment              0.114  1  0.0244 0.8790
## Mother_CEWL           14.925  1  3.2007 0.1039
## Tb_CEWL_c_mean         1.529  1  0.3279 0.5795
## msmt_temp_C_mean       5.518  1  1.1832 0.3022
## msmt_RH_percent_mean   5.851  1  1.2546 0.2889
## VPD_kPa_mean           5.739  1  1.2308 0.2932
## Treatment:Mother_CEWL  3.572  1  0.7661 0.4020
## Residuals             46.632 10

The first variable to drop from the model is humidity.

DAB_CEWL_LM_2 <- lm(data = dab_clutch_mom,
                      CEWL_g_m2h_mean ~ 
                      Treatment *
                      Mother_CEWL +
                      Tb_CEWL_c_mean +
                      msmt_temp_C_mean +
                      VPD_kPa_mean)
                     
car::vif(DAB_CEWL_LM_2, type = "predictor")
##                      GVIF Df GVIF^(1/(2*Df)) Interacts With
## Treatment        2.026678  3        1.124944    Mother_CEWL
## Mother_CEWL      2.026678  3        1.124944      Treatment
## Tb_CEWL_c_mean   1.232856  1        1.110341           --  
## msmt_temp_C_mean 4.836938  1        2.199304           --  
## VPD_kPa_mean     4.450099  1        2.109526           --  
##                                                          Other Predictors
## Treatment                  Tb_CEWL_c_mean, msmt_temp_C_mean, VPD_kPa_mean
## Mother_CEWL                Tb_CEWL_c_mean, msmt_temp_C_mean, VPD_kPa_mean
## Tb_CEWL_c_mean     Treatment, Mother_CEWL, msmt_temp_C_mean, VPD_kPa_mean
## msmt_temp_C_mean     Treatment, Mother_CEWL, Tb_CEWL_c_mean, VPD_kPa_mean
## VPD_kPa_mean     Treatment, Mother_CEWL, Tb_CEWL_c_mean, msmt_temp_C_mean
drop1(DAB_CEWL_LM_2)
## Single term deletions
## 
## Model:
## CEWL_g_m2h_mean ~ Treatment * Mother_CEWL + Tb_CEWL_c_mean + 
##     msmt_temp_C_mean + VPD_kPa_mean
##                       Df Sum of Sq    RSS    AIC
## <none>                             52.482 33.262
## Tb_CEWL_c_mean         1    7.3271 59.809 33.614
## msmt_temp_C_mean       1    0.2636 52.746 31.352
## VPD_kPa_mean           1    0.0840 52.566 31.291
## Treatment:Mother_CEWL  1    9.8849 62.367 34.368
car::Anova(DAB_CEWL_LM_2, type = "2")
## Anova Table (Type II tests)
## 
## Response: CEWL_g_m2h_mean
##                       Sum Sq Df F value  Pr(>F)  
## Treatment              7.273  1  1.5243 0.24269  
## Mother_CEWL           29.399  1  6.1619 0.03045 *
## Tb_CEWL_c_mean         7.327  1  1.5357 0.24105  
## msmt_temp_C_mean       0.264  1  0.0553 0.81847  
## VPD_kPa_mean           0.084  1  0.0176 0.89682  
## Treatment:Mother_CEWL  9.885  1  2.0718 0.17789  
## Residuals             52.482 11                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DAB_CEWL_LM_3 <- lm(data = dab_clutch_mom,
                      CEWL_g_m2h_mean ~ 
                      Treatment *
                      Mother_CEWL +
                      Tb_CEWL_c_mean +
                      msmt_temp_C_mean)
                     
drop1(DAB_CEWL_LM_3)
## Single term deletions
## 
## Model:
## CEWL_g_m2h_mean ~ Treatment * Mother_CEWL + Tb_CEWL_c_mean + 
##     msmt_temp_C_mean
##                       Df Sum of Sq    RSS    AIC
## <none>                             52.566 31.291
## Tb_CEWL_c_mean         1    8.0957 60.662 31.869
## msmt_temp_C_mean       1    0.2507 52.817 29.376
## Treatment:Mother_CEWL  1   10.4629 63.029 32.558
car::Anova(DAB_CEWL_LM_3, type = "2")
## Anova Table (Type II tests)
## 
## Response: CEWL_g_m2h_mean
##                       Sum Sq Df F value  Pr(>F)  
## Treatment              7.365  1  1.6812 0.21914  
## Mother_CEWL           35.111  1  8.0153 0.01514 *
## Tb_CEWL_c_mean         8.096  1  1.8481 0.19900  
## msmt_temp_C_mean       0.251  1  0.0572 0.81496  
## Treatment:Mother_CEWL 10.463  1  2.3885 0.14818  
## Residuals             52.566 12                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DAB_CEWL_LM_4 <- lm(data = dab_clutch_mom,
                      CEWL_g_m2h_mean ~ 
                      Treatment *
                      Mother_CEWL +
                      Tb_CEWL_c_mean)
                     
drop1(DAB_CEWL_LM_4)
## Single term deletions
## 
## Model:
## CEWL_g_m2h_mean ~ Treatment * Mother_CEWL + Tb_CEWL_c_mean
##                       Df Sum of Sq    RSS    AIC
## <none>                             52.817 29.376
## Tb_CEWL_c_mean         1    10.033 62.850 30.507
## Treatment:Mother_CEWL  1    11.926 64.743 31.041
car::Anova(DAB_CEWL_LM_4, type = "2")
## Anova Table (Type II tests)
## 
## Response: CEWL_g_m2h_mean
##                       Sum Sq Df F value   Pr(>F)   
## Treatment              9.769  1  2.4045 0.144985   
## Mother_CEWL           41.815  1 10.2921 0.006859 **
## Tb_CEWL_c_mean        10.033  1  2.4694 0.140096   
## Treatment:Mother_CEWL 11.926  1  2.9355 0.110382   
## Residuals             52.817 13                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DAB_CEWL_LM_5 <- lm(data = dab_clutch_mom,
                      CEWL_g_m2h_mean ~ 
                      Treatment +
                      Mother_CEWL +
                      Tb_CEWL_c_mean)
                     
drop1(DAB_CEWL_LM_5)
## Single term deletions
## 
## Model:
## CEWL_g_m2h_mean ~ Treatment + Mother_CEWL + Tb_CEWL_c_mean
##                Df Sum of Sq     RSS    AIC
## <none>                       64.743 31.041
## Treatment       1     9.769  74.512 31.571
## Mother_CEWL     1    41.815 106.559 38.010
## Tb_CEWL_c_mean  1    14.512  79.256 32.682
car::Anova(DAB_CEWL_LM_5, type = "2")
## Anova Table (Type II tests)
## 
## Response: CEWL_g_m2h_mean
##                Sum Sq Df F value  Pr(>F)   
## Treatment       9.769  1  2.1124 0.16815   
## Mother_CEWL    41.815  1  9.0421 0.00942 **
## Tb_CEWL_c_mean 14.512  1  3.1381 0.09824 . 
## Residuals      64.743 14                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DAB_CEWL_LM_6 <- lm(data = dab_clutch_mom,
                      CEWL_g_m2h_mean ~ 
                      Mother_CEWL +
                      Tb_CEWL_c_mean)

drop1(DAB_CEWL_LM_6)
## Single term deletions
## 
## Model:
## CEWL_g_m2h_mean ~ Mother_CEWL + Tb_CEWL_c_mean
##                Df Sum of Sq     RSS    AIC
## <none>                       74.512 31.571
## Mother_CEWL     1    43.969 118.481 37.919
## Tb_CEWL_c_mean  1    15.255  89.767 32.923
car::Anova(DAB_CEWL_LM_6, type = "2")
## Anova Table (Type II tests)
## 
## Response: CEWL_g_m2h_mean
##                Sum Sq Df F value   Pr(>F)   
## Mother_CEWL    43.969  1  8.8513 0.009439 **
## Tb_CEWL_c_mean 15.255  1  3.0709 0.100117   
## Residuals      74.512 15                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DAB_CEWL_LM_7 <- lm(data = dab_clutch_mom,
                      CEWL_g_m2h_mean ~ 
                      Mother_CEWL)

car::Anova(DAB_CEWL_LM_7, type = "2")
## Anova Table (Type II tests)
## 
## Response: CEWL_g_m2h_mean
##             Sum Sq Df F value   Pr(>F)   
## Mother_CEWL 52.881  1  9.4255 0.007324 **
## Residuals   89.767 16                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Null model

DAB_CEWL_LM_null <- lm(data = dab_clutch_mom,
                      CEWL_g_m2h_mean ~ 1)

Model Selection

clutch_cewl_models_all <- list(DAB_CEWL_LM_2, 
                               DAB_CEWL_LM_3, 
                               DAB_CEWL_LM_4, 
                               DAB_CEWL_LM_5, 
                               DAB_CEWL_LM_6, 
                               DAB_CEWL_LM_7, 
                               DAB_CEWL_LM_null)

#specify model names
clutch_cewl_mod_names_all <- c('model2', 'model3',
                               'model4', 'model5',
                               'model6', 'model7',
                               'modelnull')

#calculate AIC of each model
clutch_cewl_AICc_all <- data.frame(aictab(cand.set = clutch_cewl_models_all, 
                                 modnames = clutch_cewl_mod_names_all))
clutch_cewl_AICc_all
##    Modnames K      AICc  Delta_AICc     ModelLik       AICcWt        LL
## 6    model7 3  87.71926  0.00000000 1.0000000000 0.3469297055 -40.00249
## 5    model6 4  87.72938  0.01011992 0.9949528202 0.3451786889 -38.32623
## 4    model5 5  89.12286  1.40359817 0.4956927089 0.1719705255 -37.06143
## 3    model4 6  90.09444  2.37517755 0.3049556953 0.1057981896 -35.22904
## 7 modelnull 2  93.14190  5.42264649 0.0664488206 0.0230530698 -44.17095
## 2    model3 7  95.57243  7.85316762 0.0197108938 0.0068382946 -35.18621
## 1    model2 8 102.34363 14.62436885 0.0006673577 0.0002315262 -35.17181
##      Cum.Wt
## 6 0.3469297
## 5 0.6921084
## 4 0.8640789
## 3 0.9698771
## 7 0.9929302
## 2 0.9997685
## 1 1.0000000

So the top three equivalent models are 7, 6, & 5. 6/7 are almost exactly the same, with 5 trailing a bit behind.

Model 7 is clutch avg CEWL ~ mother CEWL. Model 6 adds clutch average body temperature at the time of measurement. Model 5 adds Tb and mother’s treatment. However, Tb and mother’s treatment did not have significant effects.

Final Model

DAB_CEWL_LM_BEST <- lm(data = dab_clutch_mom,
                      CEWL_g_m2h_mean ~ 
                      Mother_CEWL)

car::Anova(DAB_CEWL_LM_BEST, type = "2")
## Anova Table (Type II tests)
## 
## Response: CEWL_g_m2h_mean
##             Sum Sq Df F value   Pr(>F)   
## Mother_CEWL 52.881  1  9.4255 0.007324 **
## Residuals   89.767 16                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(DAB_CEWL_LM_BEST)
## 
## Call:
## lm(formula = CEWL_g_m2h_mean ~ Mother_CEWL, data = dab_clutch_mom)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8832 -1.8846  0.0906  1.6421  4.8088 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   4.4366     2.5501    1.74  0.10109   
## Mother_CEWL   0.3975     0.1295    3.07  0.00732 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.369 on 16 degrees of freedom
## Multiple R-squared:  0.3707, Adjusted R-squared:  0.3314 
## F-statistic: 9.425 on 1 and 16 DF,  p-value: 0.007324

Model Assumptions

Check linear regression assumptions/conditions:

plot(DAB_CEWL_LM_BEST)

shapiro.test(residuals(DAB_CEWL_LM_BEST))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(DAB_CEWL_LM_BEST)
## W = 0.97, p-value = 0.7977

Check Tb

Tb_Ta_lm <- lm(data = dab_all, Tb_CEWL_c ~ msmt_temp_C)
summary(Tb_Ta_lm)
## 
## Call:
## lm(formula = Tb_CEWL_c ~ msmt_temp_C, data = dab_all)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0528 -0.7894 -0.1074  0.7566  2.9472 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  18.8366     3.9149   4.811 7.18e-06 ***
## msmt_temp_C   0.3643     0.1518   2.400   0.0188 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.21 on 78 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.06876,    Adjusted R-squared:  0.05682 
## F-statistic: 5.759 on 1 and 78 DF,  p-value: 0.01879
DAB_Tb <- lmerTest::lmer(data = dab_all,
                      Tb_CEWL_c ~ 
                      Treatment +
                        (1|Mother_ID))
                     
anova(DAB_Tb, test = "F", type = "2")
## Type II Analysis of Variance Table with Satterthwaite's method
##             Sum Sq  Mean Sq NumDF DenDF F value Pr(>F)
## Treatment 0.045615 0.045615     1 16.51  0.0418 0.8406
emmeans(DAB_Tb, pairwise ~ Treatment)
## $emmeans
##  Treatment emmean    SE   df lower.CL upper.CL
##  C           28.3 0.312 16.3     27.6     28.9
##  W           28.2 0.276 15.7     27.6     28.8
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate    SE df t.ratio p.value
##  C - W      0.0851 0.416 16   0.204  0.8407
## 
## Degrees-of-freedom method: kenward-roger

CEWL ~ osml

CEWL_osml_lm <- lm(data = dab_all, CEWL_g_m2h ~ Plasma_Osmol_Rep_mean)
summary(CEWL_osml_lm)
## 
## Call:
## lm(formula = CEWL_g_m2h ~ Plasma_Osmol_Rep_mean, data = dab_all)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.1094 -2.1594 -0.2557  1.6502 12.0127 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)  
## (Intercept)           -3.86750    8.20803  -0.471   0.6388  
## Plasma_Osmol_Rep_mean  0.05364    0.02772   1.935   0.0566 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.814 on 78 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.04581,    Adjusted R-squared:  0.03358 
## F-statistic: 3.745 on 1 and 78 DF,  p-value: 0.05659
CEWL_osml_lm2 <- lm(data = ps_all, CEWL_g_m2h ~ Plasma_Osmol_Rep_mean)
summary(CEWL_osml_lm2)
## 
## Call:
## lm(formula = CEWL_g_m2h ~ Plasma_Osmol_Rep_mean, data = ps_all)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.0605 -2.0276 -0.0202  2.8676  5.3066 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)  
## (Intercept)           10.55677    4.82017   2.190   0.0357 *
## Plasma_Osmol_Rep_mean  0.01723    0.01514   1.138   0.2632  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.956 on 33 degrees of freedom
## Multiple R-squared:  0.03778,    Adjusted R-squared:  0.008624 
## F-statistic: 1.296 on 1 and 33 DF,  p-value: 0.2632

DAB CEWL DIFFS

DAB_CEWL_LMM_0 <- lme4::lmer(data = dab_all,
                      CEWL_g_m2h ~ 
                      Treatment *
                      Mother_CEWL +
                      Tb_CEWL_c +
                      msmt_temp_C +
                      msmt_RH_percent +
                      VPD_kPa +
                        (1|Mother_ID))
                     
car::vif(DAB_CEWL_LMM_0)
##             Treatment           Mother_CEWL             Tb_CEWL_c 
##             42.504746              5.885250              1.115189 
##           msmt_temp_C       msmt_RH_percent               VPD_kPa 
##            254.801919            288.145293            691.560970 
## Treatment:Mother_CEWL 
##             41.228742
drop1(DAB_CEWL_LMM_0)
## Single term deletions
## 
## Model:
## CEWL_g_m2h ~ Treatment * Mother_CEWL + Tb_CEWL_c + msmt_temp_C + 
##     msmt_RH_percent + VPD_kPa + (1 | Mother_ID)
##                       npar    AIC
## <none>                     416.55
## Tb_CEWL_c                1 416.73
## msmt_temp_C              1 417.66
## msmt_RH_percent          1 418.07
## VPD_kPa                  1 417.94
## Treatment:Mother_CEWL    1 415.76
anova(DAB_CEWL_LMM_0, test = "F", type = "2")
## Analysis of Variance Table
##                       npar  Sum Sq Mean Sq F value
## Treatment                1  26.583  26.583  3.4268
## Mother_CEWL              1 102.933 102.933 13.2691
## Tb_CEWL_c                1  36.748  36.748  4.7372
## msmt_temp_C              1  12.596  12.596  1.6237
## msmt_RH_percent          1   0.320   0.320  0.0413
## VPD_kPa                  1  33.636  33.636  4.3361
## Treatment:Mother_CEWL    1   6.788   6.788  0.8750
DAB_CEWL_LMM_1 <- lme4::lmer(data = dab_all,
                      CEWL_g_m2h ~ 
                      Treatment +
                      Mother_CEWL +
                      Tb_CEWL_c +
                      msmt_temp_C +
                      msmt_RH_percent +
                      VPD_kPa +
                        (1|Mother_ID))
DAB_CEWL_LMM_11 <- lme4::lmer(data = dab_all,
                      CEWL_g_m2h ~ 
                      Treatment *
                      Mother_CEWL +
                      Tb_CEWL_c +
                      msmt_temp_C +
                      VPD_kPa +
                        (1|Mother_ID))
                     
car::vif(DAB_CEWL_LMM_1)
##       Treatment     Mother_CEWL       Tb_CEWL_c     msmt_temp_C msmt_RH_percent 
##        1.467280        1.483813        1.116464      221.687480      251.824614 
##         VPD_kPa 
##      595.768696
car::vif(DAB_CEWL_LMM_11)
##             Treatment           Mother_CEWL             Tb_CEWL_c 
##             33.946352              5.530680              1.052778 
##           msmt_temp_C               VPD_kPa Treatment:Mother_CEWL 
##              3.798479              3.643235             35.785122
drop1(DAB_CEWL_LMM_11)
## Single term deletions
## 
## Model:
## CEWL_g_m2h ~ Treatment * Mother_CEWL + Tb_CEWL_c + msmt_temp_C + 
##     VPD_kPa + (1 | Mother_ID)
##                       npar    AIC
## <none>                     418.07
## Tb_CEWL_c                1 419.69
## msmt_temp_C              1 416.81
## VPD_kPa                  1 416.29
## Treatment:Mother_CEWL    1 419.24

start good models

(good = no collinearity)

DAB_CEWL_LMM_1a <- lme4::lmer(data = dab_all,
                      CEWL_g_m2h ~ 
                      Treatment + Mother_CEWL +
                      Tb_CEWL_c +
                      msmt_temp_C + # remove RH to fix multicollinearity
                      VPD_kPa +
                        (1|Mother_ID))
                     
car::vif(DAB_CEWL_LMM_1a)
##   Treatment Mother_CEWL   Tb_CEWL_c msmt_temp_C     VPD_kPa 
##    1.047027    1.409371    1.044307    3.726919    3.355360
drop1(DAB_CEWL_LMM_1a)
## Single term deletions
## 
## Model:
## CEWL_g_m2h ~ Treatment + Mother_CEWL + Tb_CEWL_c + msmt_temp_C + 
##     VPD_kPa + (1 | Mother_ID)
##             npar    AIC
## <none>           419.24
## Treatment      1 419.01
## Mother_CEWL    1 423.16
## Tb_CEWL_c      1 420.97
## msmt_temp_C    1 417.77
## VPD_kPa        1 417.24
anova(DAB_CEWL_LMM_1a, test = "F", type = "2")
## Analysis of Variance Table
##             npar Sum Sq Mean Sq F value
## Treatment      1 21.414  21.414  2.7518
## Mother_CEWL    1 82.634  82.634 10.6187
## Tb_CEWL_c      1 32.242  32.242  4.1432
## msmt_temp_C    1 11.823  11.823  1.5193
## VPD_kPa        1  0.089   0.089  0.0115
DAB_CEWL_LMM_2 <- lme4::lmer(data = dab_all,
                      CEWL_g_m2h ~ 
                        Treatment + Mother_CEWL +
                      Tb_CEWL_c +
                      msmt_temp_C +
                        (1|Mother_ID))

DAB_CEWL_LMM_2v <- lme4::lmer(data = dab_all,
                      CEWL_g_m2h ~ 
                        Treatment + Mother_CEWL +
                      Tb_CEWL_c +
                      VPD_kPa +
                        (1|Mother_ID))
                     
summary(DAB_CEWL_LMM_2v)
## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ Treatment + Mother_CEWL + Tb_CEWL_c + VPD_kPa +  
##     (1 | Mother_ID)
##    Data: dab_all
## 
## REML criterion at convergence: 400.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3257 -0.5481 -0.1215  0.3664  3.1862 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Mother_ID (Intercept) 2.889    1.700   
##  Residual              7.857    2.803   
## Number of obs: 80, groups:  Mother_ID, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) -17.5678    10.3296  -1.701
## TreatmentW   -1.3391     1.0469  -1.279
## Mother_CEWL   0.3917     0.1192   3.285
## Tb_CEWL_c     0.5949     0.2987   1.991
## VPD_kPa       3.0621     3.1028   0.987
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmnW M_CEWL T_CEWL
## TreatmentW  -0.181                     
## Mother_CEWL -0.207  0.075              
## Tb_CEWL_c   -0.749  0.001 -0.088       
## VPD_kPa     -0.556  0.179  0.090 -0.081
drop1(DAB_CEWL_LMM_2)
## Single term deletions
## 
## Model:
## CEWL_g_m2h ~ Treatment + Mother_CEWL + Tb_CEWL_c + msmt_temp_C + 
##     (1 | Mother_ID)
##             npar    AIC
## <none>           417.24
## Treatment      1 417.01
## Mother_CEWL    1 422.89
## Tb_CEWL_c      1 419.03
## msmt_temp_C    1 416.92
anova(DAB_CEWL_LMM_2, test = "F", type = "2")
## Analysis of Variance Table
##             npar Sum Sq Mean Sq F value
## Treatment      1 22.616  22.616  2.9173
## Mother_CEWL    1 87.344  87.344 11.2667
## Tb_CEWL_c      1 33.280  33.280  4.2929
## msmt_temp_C    1 12.006  12.006  1.5487
DAB_CEWL_LMM_3 <- lme4::lmer(data = dab_all,
                      CEWL_g_m2h ~ 
                        Treatment + Mother_CEWL +
                      Tb_CEWL_c +
                        (1|Mother_ID))

DAB_CEWL_LMM_3int <- lme4::lmer(data = dab_all,
                      CEWL_g_m2h ~ 
                        Treatment * Mother_CEWL +
                      Tb_CEWL_c +
                        (1|Mother_ID))
                     
summary(DAB_CEWL_LMM_3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ Treatment + Mother_CEWL + Tb_CEWL_c + (1 | Mother_ID)
##    Data: dab_all
## 
## REML criterion at convergence: 405.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3067 -0.5956 -0.1142  0.3314  3.2341 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Mother_ID (Intercept) 2.846    1.687   
##  Residual              7.870    2.805   
## Number of obs: 80, groups:  Mother_ID, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) -11.9346     8.5779  -1.391
## TreatmentW   -1.5239     1.0257  -1.486
## Mother_CEWL   0.3812     0.1182   3.224
## Tb_CEWL_c     0.6201     0.2977   2.083
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmnW M_CEWL
## TreatmentW  -0.099              
## Mother_CEWL -0.189  0.060       
## Tb_CEWL_c   -0.959  0.016 -0.082
drop1(DAB_CEWL_LMM_3)
## Single term deletions
## 
## Model:
## CEWL_g_m2h ~ Treatment + Mother_CEWL + Tb_CEWL_c + (1 | Mother_ID)
##             npar    AIC
## <none>           416.92
## Treatment      1 417.44
## Mother_CEWL    1 424.45
## Tb_CEWL_c      1 419.59
anova(DAB_CEWL_LMM_3, test = "F", type = "2")
## Analysis of Variance Table
##             npar Sum Sq Mean Sq F value
## Treatment      1 23.621  23.621  3.0013
## Mother_CEWL    1 91.286  91.286 11.5989
## Tb_CEWL_c      1 34.153  34.153  4.3395
DAB_CEWL_LMM_4 <- lme4::lmer(data = dab_all,
                      CEWL_g_m2h ~ 
                        Mother_CEWL +
                      Tb_CEWL_c +
                        (1|Mother_ID))
                     
summary(DAB_CEWL_LMM_4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ Mother_CEWL + Tb_CEWL_c + (1 | Mother_ID)
##    Data: dab_all
## 
## REML criterion at convergence: 410
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4120 -0.5772 -0.1474  0.4340  3.2200 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Mother_ID (Intercept) 3.205    1.790   
##  Residual              7.869    2.805   
## Number of obs: 80, groups:  Mother_ID, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) -12.9081     8.6252  -1.497
## Mother_CEWL   0.3913     0.1225   3.194
## Tb_CEWL_c     0.6173     0.3000   2.058
## 
## Correlation of Fixed Effects:
##             (Intr) M_CEWL
## Mother_CEWL -0.195       
## Tb_CEWL_c   -0.960 -0.080
drop1(DAB_CEWL_LMM_4)
## Single term deletions
## 
## Model:
## CEWL_g_m2h ~ Mother_CEWL + Tb_CEWL_c + (1 | Mother_ID)
##             npar    AIC
## <none>           417.44
## Mother_CEWL    1 424.34
## Tb_CEWL_c      1 419.81
anova(DAB_CEWL_LMM_4, test = "F", type = "2")
## Analysis of Variance Table
##             npar Sum Sq Mean Sq F value
## Mother_CEWL    1 89.300  89.300 11.3486
## Tb_CEWL_c      1 33.325  33.325  4.2351
DAB_CEWL_LMM_5 <- lme4::lmer(data = dab_all,
                      CEWL_g_m2h ~ 
                        Mother_CEWL +
                        (1|Mother_ID))
                     
summary(DAB_CEWL_LMM_5)
## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ Mother_CEWL + (1 | Mother_ID)
##    Data: dab_all
## 
## REML criterion at convergence: 413.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7448 -0.5548 -0.1630  0.4702  3.1358 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Mother_ID (Intercept) 3.920    1.980   
##  Residual              7.978    2.825   
## Number of obs: 80, groups:  Mother_ID, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   4.1727     2.5782   1.618
## Mother_CEWL   0.4097     0.1308   3.131
## 
## Correlation of Fixed Effects:
##             (Intr)
## Mother_CEWL -0.976
anova(DAB_CEWL_LMM_5, test = "F", type = "2")
## Analysis of Variance Table
##             npar Sum Sq Mean Sq F value
## Mother_CEWL    1 78.209  78.209  9.8027
DAB_CEWL_LMM_NULL <- lme4::lmer(data = dab_all,
                      CEWL_g_m2h ~ 1 + (1|Mother_ID))

Model Selection

cewl_models_all <- list(DAB_CEWL_LMM_1a, DAB_CEWL_LMM_2, DAB_CEWL_LMM_2v,
                        DAB_CEWL_LMM_3, DAB_CEWL_LMM_4, DAB_CEWL_LMM_5,
                        DAB_CEWL_LMM_NULL)

#specify model names
cewl_mod_names_all <- c('model1', 
                    'model2', 'model2 VPD', 
                    'model3',
                    'model4',
                    'model5',
                    'modelnull')

#calculate AIC of each model
cewl_AICc_all <- data.frame(aictab(cand.set = cewl_models_all, 
                                 modnames = cewl_mod_names_all))
cewl_AICc_all
##     Modnames K     AICc Delta_AICc    ModelLik      AICcWt    Res.LL    Cum.Wt
## 1     model1 8 416.2621   0.000000 1.000000000 0.385961219 -199.1170 0.3859612
## 3 model2 VPD 7 416.4001   0.137955 0.933347685 0.360236010 -200.4223 0.7461972
## 4     model3 6 419.0710   2.808868 0.245505920 0.094755764 -202.9602 0.8409530
## 2     model2 7 419.1056   2.843511 0.241290072 0.093128610 -201.7750 0.9340816
## 5     model4 5 420.7823   4.520133 0.104343549 0.040272563 -204.9857 0.9743542
## 6     model5 4 422.0008   5.738702 0.056735748 0.021897799 -206.7338 0.9962520
## 7  modelnull 3 425.5311   9.269010 0.009710911 0.003748035 -209.6077 1.0000000

Okay! after some additional investigations, it seems like the model needs to include VPD.

DAB_CEWL_LMM_2v <- lme4::lmer(data = dab_all, CEWL_g_m2h ~ Treatment + Mother_CEWL + Tb_CEWL_c + VPD_kPa + (1|Mother_ID))

Final Model + emmeans

DAB_CEWL_LMM_BEST_simple <- lmerTest::lmer(data = dab_all,
                      CEWL_g_m2h ~ 
                        Treatment +
                        Mother_CEWL +
                        Tb_CEWL_c +
                        VPD_kPa +
                          (1|Mother_ID))

summary(DAB_CEWL_LMM_BEST_simple)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h ~ Treatment + Mother_CEWL + Tb_CEWL_c + VPD_kPa +  
##     (1 | Mother_ID)
##    Data: dab_all
## 
## REML criterion at convergence: 400.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3257 -0.5481 -0.1215  0.3664  3.1862 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Mother_ID (Intercept) 2.889    1.700   
##  Residual              7.857    2.803   
## Number of obs: 80, groups:  Mother_ID, 18
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)   
## (Intercept) -17.5678    10.3296  43.4663  -1.701  0.09614 . 
## TreatmentW   -1.3391     1.0469  14.2913  -1.279  0.22125   
## Mother_CEWL   0.3917     0.1192  14.4439   3.285  0.00523 **
## Tb_CEWL_c     0.5949     0.2987  74.7079   1.991  0.05009 . 
## VPD_kPa       3.0621     3.1028  21.2580   0.987  0.33480   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmnW M_CEWL T_CEWL
## TreatmentW  -0.181                     
## Mother_CEWL -0.207  0.075              
## Tb_CEWL_c   -0.749  0.001 -0.088       
## VPD_kPa     -0.556  0.179  0.090 -0.081
anova(DAB_CEWL_LMM_BEST_simple, test = "F", type = "2")
## Type II Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## Treatment   12.855  12.855     1 14.291  1.6361 0.221245   
## Mother_CEWL 84.809  84.809     1 14.444 10.7935 0.005227 **
## Tb_CEWL_c   31.160  31.160     1 74.708  3.9657 0.050093 . 
## VPD_kPa      7.653   7.653     1 21.258  0.9739 0.334802   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(DAB_CEWL_LMM_BEST_simple)
##            R2m       R2c
## [1,] 0.3075859 0.4937201
emmeans(DAB_CEWL_LMM_BEST_simple, pairwise ~ Treatment)
## $emmeans
##  Treatment emmean    SE   df lower.CL upper.CL
##  C           12.8 0.784 13.9    11.10     14.5
##  W           11.4 0.683 13.3     9.97     12.9
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate   SE   df t.ratio p.value
##  C - W        1.34 1.05 13.7   1.277  0.2228
## 
## Degrees-of-freedom method: kenward-roger
DAB_CEWL_emmeans <- emmeans(DAB_CEWL_LMM_BEST_simple, "Treatment")
DAB_CEWL_emmeans_CI <- tidy(confint(DAB_CEWL_emmeans))
DAB_CEWL_LMM_BEST_full <- lmerTest::lmer(data = dab_all,
                      CEWL_g_m2h ~ 
                      Treatment +
                      Mother_CEWL +
                      Tb_CEWL_c +
                      msmt_temp_C +
                      VPD_kPa +
                        (1|Mother_ID))

summary(DAB_CEWL_LMM_BEST_full)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h ~ Treatment + Mother_CEWL + Tb_CEWL_c + msmt_temp_C +  
##     VPD_kPa + (1 | Mother_ID)
##    Data: dab_all
## 
## REML criterion at convergence: 398.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3517 -0.5491 -0.1371  0.3740  3.1195 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Mother_ID (Intercept) 3.263    1.806   
##  Residual              7.782    2.790   
## Number of obs: 80, groups:  Mother_ID, 18
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept) -30.0691    19.4637  21.4560  -1.545   0.1370  
## TreatmentW   -1.2540     1.0912  13.4906  -1.149   0.2704  
## Mother_CEWL   0.3317     0.1455  16.1331   2.280   0.0366 *
## Tb_CEWL_c     0.5450     0.3045  73.9954   1.790   0.0776 .
## msmt_temp_C   0.8627     1.1108  23.6049   0.777   0.4451  
## VPD_kPa      -0.6125     5.7224  26.0644  -0.107   0.9156  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmnW M_CEWL T_CEWL msm__C
## TreatmentW  -0.180                            
## Mother_CEWL  0.345  0.012                     
## Tb_CEWL_c   -0.247 -0.016  0.022              
## msmt_temp_C -0.842  0.097 -0.527 -0.175       
## VPD_kPa      0.528  0.019  0.480  0.103 -0.830
r.squaredGLMM(DAB_CEWL_LMM_BEST_full)
##            R2m      R2c
## [1,] 0.3118644 0.515141
anova(DAB_CEWL_LMM_BEST_full, test = "F", type = "2")
## Type II Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Treatment   10.277  10.277     1 13.491  1.3207 0.27044  
## Mother_CEWL 40.439  40.439     1 16.133  5.1965 0.03656 *
## Tb_CEWL_c   24.929  24.929     1 73.995  3.2034 0.07758 .
## msmt_temp_C  4.693   4.693     1 23.605  0.6031 0.44512  
## VPD_kPa      0.089   0.089     1 26.064  0.0115 0.91558  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model Assumptions

plot(DAB_CEWL_LMM_BEST_simple)

shapiro.test(residuals(DAB_CEWL_LMM_BEST_simple))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(DAB_CEWL_LMM_BEST_simple)
## W = 0.8914, p-value = 5.377e-06

looks fine, but not great

Partial Regression

plot(effect("Mother_CEWL", DAB_CEWL_LMM_BEST_simple))

plot(effect("Tb_CEWL_c", DAB_CEWL_LMM_BEST_simple))

plot(effect("VPD_kPa", DAB_CEWL_LMM_BEST_simple))

plot(effect("Treatment", DAB_CEWL_LMM_BEST_simple))

mother_effects_data <- summary(effect("Mother_CEWL",
                              DAB_CEWL_LMM_BEST_simple))
mother_effects_df <- data_frame(
  meanval = mother_effects_data$effect,
  lowCI = mother_effects_data$lower,
  highCI = mother_effects_data$upper
) %>% 
  mutate(Mother_CEWL = c(10, 15, 19, 23, 28))

mother_effects_df
## # A tibble: 5 × 4
##     meanval     lowCI    highCI Mother_CEWL
##   <dbl[1d]> <dbl[1d]> <dbl[1d]>       <dbl>
## 1      8.39      5.97      10.8          10
## 2     10.3       8.92      11.8          15
## 3     11.9      10.9       12.9          19
## 4     13.5      12.1       14.8          23
## 5     15.4      13.1       17.8          28
Tb_effects_data <- summary(effect("Tb_CEWL_c",
                        DAB_CEWL_LMM_BEST_simple))
Tb_effects_df <- data_frame(
  meanval = Tb_effects_data$effect,
  lowCI = Tb_effects_data$lower,
  highCI = Tb_effects_data$upper
) %>% 
  mutate(Tb_CEWL_c = c(26, 27, 29, 30, 31))

Tb_effects_df
## # A tibble: 5 × 4
##     meanval     lowCI    highCI Tb_CEWL_c
##   <dbl[1d]> <dbl[1d]> <dbl[1d]>     <dbl>
## 1      10.7      9.02      12.4        26
## 2      11.3     10.0       12.5        27
## 3      12.5     11.4       13.6        29
## 4      13.1     11.6       14.5        30
## 5      13.7     11.7       15.6        31
VPD_effects_data <- summary(effect("VPD_kPa",
                        DAB_CEWL_LMM_BEST_simple))
VPD_effects_df <- data_frame(
  meanval = VPD_effects_data$effect,
  lowCI = VPD_effects_data$lower,
  highCI = VPD_effects_data$upper
) %>% 
  mutate(VPD_kPa = c(1.6, 1.8, 1.9, 2.1, 2.3))

VPD_effects_df
## # A tibble: 5 × 4
##     meanval     lowCI    highCI VPD_kPa
##   <dbl[1d]> <dbl[1d]> <dbl[1d]>   <dbl>
## 1      10.9      8.38      13.4     1.6
## 2      11.5     10.0       13.0     1.8
## 3      11.8     10.7       12.9     1.9
## 4      12.4     11.1       13.7     2.1
## 5      13.0     10.8       15.3     2.3

Check PS

PS_CEWL_LMM_mother_simp <- lmerTest::lmer(data = ps_all_mother,
                      CEWL_g_m2h ~ 
                        Mother_CEWL +
                          (1|Mother_ID))


summary(PS_CEWL_LMM_mother_simp)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h ~ Mother_CEWL + (1 | Mother_ID)
##    Data: ps_all_mother
## 
## REML criterion at convergence: 175.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.75703 -0.72379  0.01135  0.73227  1.97394 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Mother_ID (Intercept) 0.000    0.000   
##  Residual              9.024    3.004   
## Number of obs: 35, groups:  Mother_ID, 8
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) 17.34677    2.91566 33.00000   5.950 1.12e-06 ***
## Mother_CEWL -0.06766    0.14575 33.00000  -0.464    0.646    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Mother_CEWL -0.985
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
anova(PS_CEWL_LMM_mother_simp, test = "F", type = "2")
## Type II Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Mother_CEWL 1.9443  1.9443     1    33  0.2155 0.6456
PS_CEWL_LMM_mother <- lmerTest::lmer(data = ps_all_mother,
                      CEWL_g_m2h ~ 
                        Treatment +
                        Mother_CEWL +
                        Tb_CEWL_c +
                        VPD_kPa +
                          (1|Mother_ID))


summary(PS_CEWL_LMM_mother)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h ~ Treatment + Mother_CEWL + Tb_CEWL_c + VPD_kPa +  
##     (1 | Mother_ID)
##    Data: ps_all_mother
## 
## REML criterion at convergence: 166.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.52396 -0.68095 -0.05939  0.69543  1.77997 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Mother_ID (Intercept) 0.000    0.000   
##  Residual              9.566    3.093   
## Number of obs: 35, groups:  Mother_ID, 8
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)
## (Intercept) 15.22859   21.21809 30.00000   0.718    0.478
## TreatmentW  -1.14817    1.19906 30.00000  -0.958    0.346
## Mother_CEWL -0.10936    0.16498 30.00000  -0.663    0.512
## Tb_CEWL_c   -0.01035    0.53759 30.00000  -0.019    0.985
## VPD_kPa      1.87699    6.70663 30.00000   0.280    0.781
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmnW M_CEWL T_CEWL
## TreatmentW  -0.339                     
## Mother_CEWL -0.293  0.279              
## Tb_CEWL_c   -0.722  0.374 -0.046       
## VPD_kPa     -0.667 -0.001  0.263 -0.010
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
anova(PS_CEWL_LMM_mother, test = "F", type = "2")
## Type II Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treatment   8.7715  8.7715     1    30  0.9169 0.3459
## Mother_CEWL 4.2030  4.2030     1    30  0.4394 0.5125
## Tb_CEWL_c   0.0035  0.0035     1    30  0.0004 0.9848
## VPD_kPa     0.7493  0.7493     1    30  0.0783 0.7815
r.squaredGLMM(PS_CEWL_LMM_mother)
##             R2m        R2c
## [1,] 0.03773164 0.03773164
emmeans(PS_CEWL_LMM_mother, pairwise ~ Treatment)
## $emmeans
##  Treatment emmean    SE   df lower.CL upper.CL
##  C           16.5 0.770 3.10     14.1     18.9
##  W           15.4 0.865 3.36     12.8     18.0
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate   SE   df t.ratio p.value
##  C - W        1.15 1.24 3.26   0.924  0.4189
## 
## Degrees-of-freedom method: kenward-roger
PS_CEWL_LMM_mother2 <- lme4::lmer(data = ps_all_mother,
                      CEWL_g_m2h ~ 
                        Treatment +
                        Mother_CEWL +
                        Tb_CEWL_c +
                        VPD_kPa +
                          (1|Mother_ID))
drop1(PS_CEWL_LMM_mother2)
## Single term deletions
## 
## Model:
## CEWL_g_m2h ~ Treatment + Mother_CEWL + Tb_CEWL_c + VPD_kPa + 
##     (1 | Mother_ID)
##             npar    AIC
## <none>           186.97
## Treatment      1 186.02
## Mother_CEWL    1 185.48
## Tb_CEWL_c      1 184.97
## VPD_kPa        1 185.06
PS_CEWL_LMM_mother3 <- lme4::lmer(data = ps_all_mother,
                      CEWL_g_m2h ~ 
                        Treatment +
                        Mother_CEWL +
                        VPD_kPa +
                          (1|Mother_ID))
drop1(PS_CEWL_LMM_mother3)
## Single term deletions
## 
## Model:
## CEWL_g_m2h ~ Treatment + Mother_CEWL + VPD_kPa + (1 | Mother_ID)
##             npar    AIC
## <none>           184.97
## Treatment      1 184.17
## Mother_CEWL    1 183.48
## VPD_kPa        1 183.06
anova(PS_CEWL_LMM_mother3)
## Analysis of Variance Table
##             npar Sum Sq Mean Sq F value
## Treatment      1 6.4075  6.4075  0.6921
## Mother_CEWL    1 5.5940  5.5940  0.6043
## VPD_kPa        1 0.7484  0.7484  0.0808

PUB FIG

CEWL CLUTCH ~ MOTHER

ggplot(
  dab_clutch_mom,
  aes(
    x = Mother_CEWL,
    y = CEWL_g_m2h_mean
  )
) +
  geom_abline(
    slope = 1, intercept = 0, 
    linetype = "dashed", color = "darkgrey"
  ) +
  geom_pointrange(
    aes(ymin = CEWL_g_m2h_mean-CEWL_g_m2h_sd,
        ymax = CEWL_g_m2h_mean+CEWL_g_m2h_sd,
        color = Treatment, shape = Treatment),
    #color = "springgreen3",
    alpha = 0.6, size = 0.2, linewidth = 0.3
  ) +
  geom_point(aes(color = Treatment, shape = Treatment)) +
  geom_smooth(
    se = F,
    formula = 'y~x',
    method = "lm",
    color = "black",
    linewidth = 1
  ) +
  scale_x_continuous(
    limits = c(0, 30),
    breaks = c(0, 10, 20, 30),
    name = expression('Mother CEWL (g '*m^-2*' '*h^-1*')')
  ) +
  scale_y_continuous(
    limits = c(0, 21),
    breaks = c(0, 10, 20),
    name = expression('Neonate CEWL (g '*m^-2*' '*h^-1*')')
  ) -> neo_mom_CEWL
neo_mom_CEWL

partreg CEWL~MOM

ggplot() +
  geom_abline(
    slope = 1, intercept = 0, 
    linetype = "dashed", color = "darkgrey"
  ) +
  geom_pointrange(
    data = dab_clutch_mom,
  aes(
    x = Mother_CEWL,
    y = CEWL_g_m2h_mean,
    ymin = CEWL_g_m2h_mean-CEWL_g_m2h_sd,
    ymax = CEWL_g_m2h_mean+CEWL_g_m2h_sd,
    color = Treatment, shape = Treatment),
  alpha = 0.6, size = 0.2, linewidth = 0.3, show.legend = FALSE
  ) +
  geom_point(
    data = dab_clutch_mom,
  aes(
    x = Mother_CEWL,
    y = CEWL_g_m2h_mean,
    color = Treatment, shape = Treatment
  )) +
  
  geom_ribbon(
    data = mother_effects_df,
    aes(x = Mother_CEWL, ymin = lowCI, ymax = highCI),
    alpha = 0.2
  ) +
  geom_line(
    data = mother_effects_df,
    aes(x = Mother_CEWL, y = meanval),
    color = "black",
    linewidth = 1
  ) +
  
  # prettys 
  scale_color_manual(
    values = c("darkgoldenrod2", "royalblue4"), 
    labels = c("C" = "Control", "W" = "Hydration"), 
    name = NULL
  ) +
  scale_shape_manual(
    values = c(15, 17), 
    labels = c("C" = "Control", "W" = "Hydration"),
    name = NULL
  ) +
  theme(
    legend.position = "inside",
    legend.position.inside = c(0.75, 0.15),
    legend.background = element_blank()
  ) +
  # annotate(
  #   geom = "text", label = expression(''*R^2*'=0.5'),
  #   x = 2, y = 13, size = 3
  # ) +
  scale_x_continuous(
    limits = c(0, 30),
    breaks = c(0, 10, 20, 30),
    name = expression('Mother CEWL (g '*m^-2*' '*h^-1*')')
  ) +
  scale_y_continuous(
    limits = c(0, 21),
    breaks = c(0, 10, 20),
    name = expression('Neonate CEWL (g '*m^-2*' '*h^-1*')')
  ) -> neo_mom_CEWL_partreg
neo_mom_CEWL_partreg

export it:

ggsave(
  filename = "Fig4_neonate_mother_CEWL.pdf",
  plot = neo_mom_CEWL_partreg,
  path = "./results_figures",
  device = "pdf",
  dpi = 600,
  units = "mm",
  width = 150,
  height = 100
)
ggsave(
  filename = "Fig4_neonate_mother_CEWL.tiff",
  plot = neo_mom_CEWL_partreg,
  path = "./results_figures",
  device = "tiff",
  dpi = 600,
  units = "mm",
  width = 150,
  height = 100
)

partreg CEWL~Tb

ggplot() +
  geom_pointrange(
    data = dab_clutch_mom,
  aes(
    x = Tb_CEWL_c_mean,
    y = CEWL_g_m2h_mean,
    ymin = CEWL_g_m2h_mean-CEWL_g_m2h_sd,
    ymax = CEWL_g_m2h_mean+CEWL_g_m2h_sd,
    color = Treatment, shape = Treatment),
  alpha = 0.6, size = 0.2, linewidth = 0.3
  ) +
  geom_point(
    data = dab_clutch_mom,
  aes(
    x = Tb_CEWL_c_mean,
    y = CEWL_g_m2h_mean,
    color = Treatment, shape = Treatment
  )) +
  
  geom_ribbon(
    data = Tb_effects_df,
    aes(x = Tb_CEWL_c, ymin = lowCI, ymax = highCI),
    alpha = 0.2
  ) +
  geom_line(
    data = Tb_effects_df,
    aes(x = Tb_CEWL_c, y = meanval),
    color = "black",
    linewidth = 1
  ) +
  
  # prettys 
  scale_color_manual(
    values = c("darkgoldenrod2", "royalblue4"), 
    labels = c("C" = "Control", "W" = "Hydration"), 
    name = NULL
  ) +
  scale_shape_manual(
    values = c(15, 17), 
    labels = c("C" = "Control", "W" = "Hydration"),
    name = NULL
  ) +
  theme(
    legend.position = "inside",
    legend.position.inside = c(0.1, 0.1),
    legend.background = element_blank()
  ) +
  # scale_x_continuous(
  #   limits = c(26, 30),
  #   breaks = c(0, 10, 20, 30),
  #   name = expression('Mother CEWL (g '*m^-2*' '*h^-1*')')
  # ) +
  scale_y_continuous(
    #limits = c(0, 21),
    breaks = c(0, 10, 20),
    name = expression('Neonate CEWL (g '*m^-2*' '*h^-1*')')
  ) -> neo_Tb_CEWL_partreg
neo_Tb_CEWL_partreg

partreg CEWL~VPD

ggplot() +
  geom_pointrange(
    data = dab_clutch_mom,
  aes(
    x = VPD_kPa_mean,
    y = CEWL_g_m2h_mean,
    ymin = CEWL_g_m2h_mean-CEWL_g_m2h_sd,
    ymax = CEWL_g_m2h_mean+CEWL_g_m2h_sd,
    color = Treatment, shape = Treatment),
  alpha = 0.6, size = 0.2, linewidth = 0.3
  ) +
  geom_point(
    data = dab_clutch_mom,
  aes(
    x = VPD_kPa_mean,
    y = CEWL_g_m2h_mean,
    color = Treatment, shape = Treatment
  )) +
  
  geom_ribbon(
    data = VPD_effects_df,
    aes(x = VPD_kPa, ymin = lowCI, ymax = highCI),
    alpha = 0.2
  ) +
  geom_line(
    data = VPD_effects_df,
    aes(x = VPD_kPa, y = meanval),
    color = "black",
    linewidth = 1
  ) +
  
  # prettys 
  scale_color_manual(
    values = c("darkgoldenrod2", "royalblue4"), 
    labels = c("C" = "Control", "W" = "Hydration"), 
    name = NULL
  ) +
  scale_shape_manual(
    values = c(15, 17), 
    labels = c("C" = "Control", "W" = "Hydration"),
    name = NULL
  ) +
  theme(
    legend.position = "inside",
    legend.position.inside = c(0.1, 0.1),
    legend.background = element_blank()
  ) +
  # scale_x_continuous(
  #   limits = c(0, 30),
  #   breaks = c(0, 10, 20, 30),
  #   name = expression('Mother CEWL (g '*m^-2*' '*h^-1*')')
  # ) +
  scale_y_continuous(
    #limits = c(0, 21),
    breaks = c(0, 10, 20),
    name = expression('Neonate CEWL (g '*m^-2*' '*h^-1*')')
  ) -> neo_VPD_CEWL_partreg
neo_VPD_CEWL_partreg

DAB osml clutch ~ mother

this is what was done originally, but we now use all neonate data and do not avg by clutch

DAB_osml_LM_1 <- lm(data = dab_clutch_mom,
                      Plasma_Osmol_Rep_mean_mean ~ 
                      Treatment *
                      (Mother_osml + Mother_Days_in_Treatment) +
                      Mass_g_mean)
                     
car::vif(DAB_osml_LM_1, type = "predictor")
##                                 GVIF Df GVIF^(1/(2*Df))
## Treatment                   1.316106  5        1.027848
## Mother_osml                 9.171779  3        1.446801
## Mother_Days_in_Treatment 1004.349474  3        3.164566
## Mass_g_mean                 1.316106  1        1.147217
##                                                 Interacts With
## Treatment                Mother_osml, Mother_Days_in_Treatment
## Mother_osml                                          Treatment
## Mother_Days_in_Treatment                             Treatment
## Mass_g_mean                                               --  
##                                                          Other Predictors
## Treatment                                                     Mass_g_mean
## Mother_osml                         Mother_Days_in_Treatment, Mass_g_mean
## Mother_Days_in_Treatment                         Mother_osml, Mass_g_mean
## Mass_g_mean              Treatment, Mother_osml, Mother_Days_in_Treatment
drop1(DAB_osml_LM_1)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean_mean ~ Treatment * (Mother_osml + Mother_Days_in_Treatment) + 
##     Mass_g_mean
##                                    Df Sum of Sq    RSS    AIC
## <none>                                          1168.5 89.116
## Mass_g_mean                         1     1.917 1170.4 87.145
## Treatment:Mother_osml               1     6.144 1174.6 87.210
## Treatment:Mother_Days_in_Treatment  1   154.172 1322.7 89.347
anova(DAB_osml_LM_1)
## Analysis of Variance Table
## 
## Response: Plasma_Osmol_Rep_mean_mean
##                                    Df  Sum Sq Mean Sq F value  Pr(>F)  
## Treatment                           1  811.93  811.93  7.6434 0.01840 *
## Mother_osml                         1  392.60  392.60  3.6959 0.08081 .
## Mother_Days_in_Treatment            1  133.48  133.48  1.2566 0.28617  
## Mass_g_mean                         1    1.02    1.02  0.0096 0.92361  
## Treatment:Mother_osml               1    4.04    4.04  0.0380 0.84896  
## Treatment:Mother_Days_in_Treatment  1  154.17  154.17  1.4514 0.25359  
## Residuals                          11 1168.49  106.23                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DAB_osml_LM_2a <- lm(data = dab_clutch_mom,
                      Plasma_Osmol_Rep_mean_mean ~ 
                      Treatment * Mother_Days_in_Treatment +
                      Mother_osml + 
                      Mass_g_mean)
                     
drop1(DAB_osml_LM_2a)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean_mean ~ Treatment * Mother_Days_in_Treatment + 
##     Mother_osml + Mass_g_mean
##                                    Df Sum of Sq    RSS    AIC
## <none>                                          1174.6 87.210
## Mother_osml                         1    43.603 1218.2 85.866
## Mass_g_mean                         1     5.030 1179.7 85.287
## Treatment:Mother_Days_in_Treatment  1   152.067 1326.7 87.401
car::Anova(DAB_osml_LM_2a, type = "2")
## Anova Table (Type II tests)
## 
## Response: Plasma_Osmol_Rep_mean_mean
##                                     Sum Sq Df F value  Pr(>F)  
## Treatment                           455.34  1  4.6517 0.05201 .
## Mother_Days_in_Treatment            124.13  1  1.2681 0.28215  
## Mother_osml                          43.60  1  0.4454 0.51714  
## Mass_g_mean                           5.03  1  0.0514 0.82449  
## Treatment:Mother_Days_in_Treatment  152.07  1  1.5535 0.23640  
## Residuals                          1174.63 12                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DAB_osml_LM_2b <- lm(data = dab_clutch_mom,
                      Plasma_Osmol_Rep_mean_mean ~ 
                      Treatment *
                      (Mother_osml + Mother_Days_in_Treatment))
drop1(DAB_osml_LM_2b)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean_mean ~ Treatment * (Mother_osml + Mother_Days_in_Treatment)
##                                    Df Sum of Sq    RSS    AIC
## <none>                                          1170.4 87.145
## Treatment:Mother_osml               1     9.256 1179.7 85.287
## Treatment:Mother_Days_in_Treatment  1   152.289 1322.7 87.347
car::Anova(DAB_osml_LM_2b, type = "2")
## Anova Table (Type II tests)
## 
## Response: Plasma_Osmol_Rep_mean_mean
##                                     Sum Sq Df F value  Pr(>F)  
## Treatment                           507.00  1  5.1982 0.04168 *
## Mother_osml                          58.46  1  0.5994 0.45381  
## Mother_Days_in_Treatment            118.37  1  1.2136 0.29222  
## Treatment:Mother_osml                 9.26  1  0.0949 0.76332  
## Treatment:Mother_Days_in_Treatment  152.29  1  1.5614 0.23528  
## Residuals                          1170.41 12                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DAB_osml_LM_3 <- lm(data = dab_clutch_mom,
                      Plasma_Osmol_Rep_mean_mean ~ 
                      Treatment * Mother_Days_in_Treatment +
                      Mother_osml)
                     
drop1(DAB_osml_LM_3)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean_mean ~ Treatment * Mother_Days_in_Treatment + 
##     Mother_osml
##                                    Df Sum of Sq    RSS    AIC
## <none>                                          1179.7 85.287
## Mother_osml                         1    58.457 1238.1 84.158
## Treatment:Mother_Days_in_Treatment  1   148.060 1327.7 85.415
car::Anova(DAB_osml_LM_3, type = "2")
## Anova Table (Type II tests)
## 
## Response: Plasma_Osmol_Rep_mean_mean
##                                     Sum Sq Df F value  Pr(>F)  
## Treatment                           507.00  1  5.5872 0.03434 *
## Mother_Days_in_Treatment            133.48  1  1.4710 0.24677  
## Mother_osml                          58.46  1  0.6442 0.43662  
## Treatment:Mother_Days_in_Treatment  148.06  1  1.6316 0.22382  
## Residuals                          1179.66 13                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DAB_osml_LM_4 <- lm(data = dab_clutch_mom,
                      Plasma_Osmol_Rep_mean_mean ~ 
                      Treatment * Mother_Days_in_Treatment)
                     
drop1(DAB_osml_LM_4)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean_mean ~ Treatment * Mother_Days_in_Treatment
##                                    Df Sum of Sq    RSS    AIC
## <none>                                          1238.1 84.158
## Treatment:Mother_Days_in_Treatment  1    225.53 1463.7 85.170
car::Anova(DAB_osml_LM_4, type = "2")
## Anova Table (Type II tests)
## 
## Response: Plasma_Osmol_Rep_mean_mean
##                                     Sum Sq Df F value   Pr(>F)   
## Treatment                           828.57  1  9.3690 0.008465 **
## Mother_Days_in_Treatment            390.16  1  4.4117 0.054293 . 
## Treatment:Mother_Days_in_Treatment  225.53  1  2.5502 0.132602   
## Residuals                          1238.12 14                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DAB_osml_LM_5 <- lm(data = dab_clutch_mom,
                      Plasma_Osmol_Rep_mean_mean ~ 
                      Treatment + Mother_Days_in_Treatment)
                     
drop1(DAB_osml_LM_5)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean_mean ~ Treatment + Mother_Days_in_Treatment
##                          Df Sum of Sq    RSS    AIC
## <none>                                1463.7 85.170
## Treatment                 1    828.57 2292.2 91.244
## Mother_Days_in_Treatment  1    390.16 1853.8 87.423
car::Anova(DAB_osml_LM_5, type = "2")
## Anova Table (Type II tests)
## 
## Response: Plasma_Osmol_Rep_mean_mean
##                           Sum Sq Df F value  Pr(>F)  
## Treatment                 828.57  1  8.4914 0.01069 *
## Mother_Days_in_Treatment  390.16  1  3.9985 0.06399 .
## Residuals                1463.65 15                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DAB_osml_LM_6 <- lm(data = dab_clutch_mom,
                      Plasma_Osmol_Rep_mean_mean ~ 
                      Treatment)

car::Anova(DAB_osml_LM_6, type = "2")
## Anova Table (Type II tests)
## 
## Response: Plasma_Osmol_Rep_mean_mean
##            Sum Sq Df F value  Pr(>F)  
## Treatment  811.93  1  7.0077 0.01757 *
## Residuals 1853.80 16                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DAB_osml_LM_null <- lm(data = dab_clutch_mom, Plasma_Osmol_Rep_mean_mean ~ 1)

Model Selection

clutch_osmol_models_all <- list(DAB_osml_LM_1, DAB_osml_LM_2a, DAB_osml_LM_3, DAB_osml_LM_4, DAB_osml_LM_5, DAB_osml_LM_6, DAB_osml_LM_null)

#specify model names
clutch_osmol_mod_names_all <- c('model1', 
                    'model2', 
                    'model3',
                    'model4',
                    'model5',
                    'model6',
                    'modelnull')

#calculate AIC of each model
clutch_osmol_AICc_all <- data.frame(aictab(cand.set = clutch_osmol_models_all, 
                                 modnames = clutch_osmol_mod_names_all))
clutch_osmol_AICc_all
##    Modnames K     AICc Delta_AICc    ModelLik       AICcWt        LL    Cum.Wt
## 5    model5 4 141.3284  0.0000000 1.000000000 4.028689e-01 -65.12573 0.4028689
## 6    model6 3 142.2193  0.8908989 0.640536320 2.580522e-01 -67.25250 0.6609211
## 4    model4 5 142.2394  0.9109883 0.634134522 2.554731e-01 -63.61969 0.9163942
## 7 modelnull 2 145.8433  4.5149020 0.104616812 4.214686e-02 -70.52165 0.9585410
## 3    model3 6 146.0052  4.6767763 0.096483027 3.887001e-02 -63.18440 0.9974111
## 2    model2 7 151.4919 10.1635002 0.006209033 2.501427e-03 -63.14595 0.9999125
## 1    model1 8 158.1975 16.8691074 0.000217230 8.751523e-05 -63.09875 1.0000000

The best top/equivalent models are 5 (tmt + mother days in tmt), 6 (tmt only), and 4 (tmt * mother days in tmt).

Final Model

DAB_osml_LM_BEST <- lm(data = dab_clutch_mom,
                       Plasma_Osmol_Rep_mean_mean ~ 
                          Treatment)

car::Anova(DAB_osml_LM_BEST, type = "2")
## Anova Table (Type II tests)
## 
## Response: Plasma_Osmol_Rep_mean_mean
##            Sum Sq Df F value  Pr(>F)  
## Treatment  811.93  1  7.0077 0.01757 *
## Residuals 1853.80 16                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(DAB_osml_LM_BEST)
## 
## Call:
## lm(formula = Plasma_Osmol_Rep_mean_mean ~ Treatment, data = dab_clutch_mom)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -20.3615  -8.7975   0.3219   7.8067  20.5719 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  304.095      3.806  79.906   <2e-16 ***
## TreatmentW   -13.516      5.106  -2.647   0.0176 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.76 on 16 degrees of freedom
## Multiple R-squared:  0.3046, Adjusted R-squared:  0.2611 
## F-statistic: 7.008 on 1 and 16 DF,  p-value: 0.01757
DAB_osml_LM_BEST_add <- lm(data = dab_clutch_mom,
                       Plasma_Osmol_Rep_mean_mean ~ 
                          Treatment + Mother_Days_in_Treatment)

car::Anova(DAB_osml_LM_BEST_add, type = "2")
## Anova Table (Type II tests)
## 
## Response: Plasma_Osmol_Rep_mean_mean
##                           Sum Sq Df F value  Pr(>F)  
## Treatment                 828.57  1  8.4914 0.01069 *
## Mother_Days_in_Treatment  390.16  1  3.9985 0.06399 .
## Residuals                1463.65 15                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(DAB_osml_LM_BEST_add)
## 
## Call:
## lm(formula = Plasma_Osmol_Rep_mean_mean ~ Treatment + Mother_Days_in_Treatment, 
##     data = dab_clutch_mom)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.0217  -7.4365   0.6639   5.2247  15.3004 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              295.9690     5.3582  55.236   <2e-16 ***
## TreatmentW               -13.6553     4.6861  -2.914   0.0107 *  
## Mother_Days_in_Treatment   0.9287     0.4644   2.000   0.0640 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.878 on 15 degrees of freedom
## Multiple R-squared:  0.4509, Adjusted R-squared:  0.3777 
## F-statistic:  6.16 on 2 and 15 DF,  p-value: 0.01115
DAB_osml_LM_BEST_int <- lm(data = dab_clutch_mom,
                       Plasma_Osmol_Rep_mean_mean ~ 
                          Treatment * Mother_Days_in_Treatment)

car::Anova(DAB_osml_LM_BEST_int, type = "2")
## Anova Table (Type II tests)
## 
## Response: Plasma_Osmol_Rep_mean_mean
##                                     Sum Sq Df F value   Pr(>F)   
## Treatment                           828.57  1  9.3690 0.008465 **
## Mother_Days_in_Treatment            390.16  1  4.4117 0.054293 . 
## Treatment:Mother_Days_in_Treatment  225.53  1  2.5502 0.132602   
## Residuals                          1238.12 14                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(DAB_osml_LM_BEST_int)
## 
## Call:
## lm(formula = Plasma_Osmol_Rep_mean_mean ~ Treatment * Mother_Days_in_Treatment, 
##     data = dab_clutch_mom)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.708  -8.931  -0.686   7.514  11.664 
## 
## Coefficients:
##                                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                         290.9261     5.9995  48.492   <2e-16 ***
## TreatmentW                           -0.9141     9.1412  -0.100   0.9218    
## Mother_Days_in_Treatment              1.5050     0.5707   2.637   0.0195 *  
## TreatmentW:Mother_Days_in_Treatment  -1.4413     0.9026  -1.597   0.1326    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.404 on 14 degrees of freedom
## Multiple R-squared:  0.5355, Adjusted R-squared:  0.436 
## F-statistic: 5.381 on 3 and 14 DF,  p-value: 0.01128

Model Assumptions

plot(DAB_osml_LM_BEST)

plot(DAB_osml_LM_BEST_add)

plot(DAB_osml_LM_BEST_int)

DAB OSML DIFFS

first, check whether tmt and mother osml are collinear

tmt_momosmol <- lm(data = dab_all, Mother_osml ~ Treatment)
summary(tmt_momosmol)
## 
## Call:
## lm(formula = Mother_osml ~ Treatment, data = dab_all)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.093  -6.093   0.366  11.240  39.490 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  317.760      2.632 120.719   <2e-16 ***
## TreatmentW    -8.875      3.493  -2.541    0.013 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.57 on 79 degrees of freedom
## Multiple R-squared:  0.07555,    Adjusted R-squared:  0.06385 
## F-statistic: 6.457 on 1 and 79 DF,  p-value: 0.01301

Yep, as expected, so I think I need to run one model selection with tmt and one with mother osml??

Tmt + Mom Osml

DAB_OSMOL_LMM_1 <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      Treatment * 
  (Mother_osml+Mother_Days_in_Treatment+live_neonates+total_embryos) + 
                        Mass_g +
                        (1|Mother_ID))

car::vif(DAB_OSMOL_LMM_1)                     
##                          Treatment                        Mother_osml 
##                         933.860788                           4.451064 
##           Mother_Days_in_Treatment                      live_neonates 
##                           4.156945                          21.171490 
##                      total_embryos                             Mass_g 
##                          15.974730                           1.506775 
##              Treatment:Mother_osml Treatment:Mother_Days_in_Treatment 
##                         984.238714                          17.840094 
##            Treatment:live_neonates            Treatment:total_embryos 
##                         383.663677                         290.979234
drop1(DAB_OSMOL_LMM_1)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean ~ Treatment * (Mother_osml + Mother_Days_in_Treatment + 
##     live_neonates + total_embryos) + Mass_g + (1 | Mother_ID)
##                                    npar    AIC
## <none>                                  653.29
## Mass_g                                1 651.64
## Treatment:Mother_osml                 1 652.27
## Treatment:Mother_Days_in_Treatment    1 652.13
## Treatment:live_neonates               1 654.19
## Treatment:total_embryos               1 653.14
DAB_OSMOL_LMM_1a <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      Treatment * 
                        (live_neonates + total_embryos) + 
                        Mother_osml + Mother_Days_in_Treatment +
                        Mass_g +
                        (1|Mother_ID))

car::vif(DAB_OSMOL_LMM_1a)                     
##                Treatment            live_neonates            total_embryos 
##                29.200490                19.241148                15.155303 
##              Mother_osml Mother_Days_in_Treatment                   Mass_g 
##                 2.698235                 2.927528                 1.480098 
##  Treatment:live_neonates  Treatment:total_embryos 
##               216.170484               156.894196
drop1(DAB_OSMOL_LMM_1a)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean ~ Treatment * (live_neonates + total_embryos) + 
##     Mother_osml + Mother_Days_in_Treatment + Mass_g + (1 | Mother_ID)
##                          npar    AIC
## <none>                        650.44
## Mother_osml                 1 648.47
## Mother_Days_in_Treatment    1 657.15
## Mass_g                      1 648.84
## Treatment:live_neonates     1 651.09
## Treatment:total_embryos     1 649.84
DAB_OSMOL_LMM_1b <- lme4::lmer(data = dab_all,
                    Plasma_Osmol_Rep_mean ~ 
                      Treatment + 
                      Mother_osml + Mother_Days_in_Treatment +
                      live_neonates + total_embryos + 
                      Mass_g +
                        (1|Mother_ID))

car::vif(DAB_OSMOL_LMM_1b)                     
##                Treatment              Mother_osml Mother_Days_in_Treatment 
##                 1.254321                 1.571342                 2.059891 
##            live_neonates            total_embryos                   Mass_g 
##                11.244267                12.120990                 1.313391
drop1(DAB_OSMOL_LMM_1b)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean ~ Treatment + Mother_osml + Mother_Days_in_Treatment + 
##     live_neonates + total_embryos + Mass_g + (1 | Mother_ID)
##                          npar    AIC
## <none>                        649.78
## Treatment                   1 652.99
## Mother_osml                 1 649.47
## Mother_Days_in_Treatment    1 653.29
## live_neonates               1 653.84
## total_embryos               1 652.06
## Mass_g                      1 648.20
anova(DAB_OSMOL_LMM_1b, test = "F", type = "2")
## Analysis of Variance Table
##                          npar  Sum Sq Mean Sq F value
## Treatment                   1 1156.42 1156.42  9.1689
## Mother_osml                 1  588.65  588.65  4.6672
## Mother_Days_in_Treatment    1  197.60  197.60  1.5667
## live_neonates               1  174.73  174.73  1.3854
## total_embryos               1  511.74  511.74  4.0574
## Mass_g                      1  105.88  105.88  0.8395

the two different clutch size assessments are too collinear, I need to choose one…

test_live <- lme4::lmer(data = dab_all,
                    Plasma_Osmol_Rep_mean ~ 
                      Treatment + 
                      Mother_osml +
                      live_neonates + 
                        (1|Mother_ID))
summary(test_live)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Plasma_Osmol_Rep_mean ~ Treatment + Mother_osml + live_neonates +  
##     (1 | Mother_ID)
##    Data: dab_all
## 
## REML criterion at convergence: 632.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7921 -0.5028 -0.0780  0.5456  3.2112 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Mother_ID (Intercept)  68.0     8.246  
##  Residual              128.3    11.327  
## Number of obs: 81, groups:  Mother_ID, 18
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)   217.1458    48.3527   4.491
## TreatmentW    -10.4905     4.9302  -2.128
## Mother_osml     0.2992     0.1491   2.007
## live_neonates  -1.0530     1.1398  -0.924
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmnW Mthr_s
## TreatmentW  -0.366              
## Mother_osml -0.977  0.300       
## live_neonts -0.171  0.095 -0.031
anova(test_live, test = "F", type = "2")
## Analysis of Variance Table
##               npar Sum Sq Mean Sq F value
## Treatment        1 995.25  995.25  7.7569
## Mother_osml      1 502.92  502.92  3.9197
## live_neonates    1 109.49  109.49  0.8534
test_total <- lme4::lmer(data = dab_all,
                    Plasma_Osmol_Rep_mean ~ 
                      Treatment + 
                      Mother_osml +
                      total_embryos + 
                        (1|Mother_ID))
summary(test_total)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Plasma_Osmol_Rep_mean ~ Treatment + Mother_osml + total_embryos +  
##     (1 | Mother_ID)
##    Data: dab_all
## 
## REML criterion at convergence: 633.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8000 -0.4869 -0.0634  0.5624  3.2034 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Mother_ID (Intercept)  71.81    8.474  
##  Residual              128.30   11.327  
## Number of obs: 81, groups:  Mother_ID, 18
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)    215.693     49.879   4.324
## TreatmentW     -10.629      5.109  -2.081
## Mother_osml      0.292      0.152   1.921
## total_embryos   -0.545      0.998  -0.546
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmnW Mthr_s
## TreatmentW  -0.386              
## Mother_osml -0.979  0.306       
## total_mbrys -0.228  0.201  0.038
anova(test_total, test = "F", type = "2")
## Analysis of Variance Table
##               npar Sum Sq Mean Sq F value
## Treatment        1 960.17  960.17  7.4839
## Mother_osml      1 484.42  484.42  3.7758
## total_embryos    1  38.25   38.25  0.2982

I think it makes sense to use the live neonates since mothers might have stopped devoting resources to slugs.

start good models

(good = no collinearity)

DAB_OSMOL_LMM_1c <- lme4::lmer(data = dab_all,
                    Plasma_Osmol_Rep_mean ~ 
                      Treatment + 
                      Mother_osml + Mother_Days_in_Treatment +
                      live_neonates +
                      Mass_g +
                        (1|Mother_ID))

car::vif(DAB_OSMOL_LMM_1c)                     
##                Treatment              Mother_osml Mother_Days_in_Treatment 
##                 1.217221                 1.508605                 1.361143 
##            live_neonates                   Mass_g 
##                 1.138326                 1.182775
drop1(DAB_OSMOL_LMM_1c)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean ~ Treatment + Mother_osml + Mother_Days_in_Treatment + 
##     live_neonates + Mass_g + (1 | Mother_ID)
##                          npar    AIC
## <none>                        652.06
## Treatment                   1 655.10
## Mother_osml                 1 652.50
## Mother_Days_in_Treatment    1 651.77
## live_neonates               1 652.82
## Mass_g                      1 652.22
anova(DAB_OSMOL_LMM_1c, test = "F", type = "2")
## Analysis of Variance Table
##                          npar Sum Sq Mean Sq F value
## Treatment                   1 954.34  954.34  7.6772
## Mother_osml                 1 481.36  481.36  3.8723
## Mother_Days_in_Treatment    1 161.88  161.88  1.3022
## live_neonates               1 145.43  145.43  1.1699
## Mass_g                      1 274.90  274.90  2.2114
DAB_OSMOL_LMM_2 <- lme4::lmer(data = dab_all,
                    Plasma_Osmol_Rep_mean ~ 
                      Treatment + 
                      Mother_osml +
                      live_neonates + 
                      Mass_g + 
                        (1|Mother_ID))

drop1(DAB_OSMOL_LMM_2)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean ~ Treatment + Mother_osml + live_neonates + 
##     Mass_g + (1 | Mother_ID)
##               npar    AIC
## <none>             651.77
## Treatment        1 653.56
## Mother_osml      1 655.31
## live_neonates    1 652.04
## Mass_g           1 652.42
anova(DAB_OSMOL_LMM_2, test = "F", type = "2")
## Analysis of Variance Table
##               npar Sum Sq Mean Sq F value
## Treatment        1 933.85  933.85  7.5226
## Mother_osml      1 470.58  470.58  3.7908
## live_neonates    1 102.96  102.96  0.8294
## Mass_g           1 326.01  326.01  2.6262
DAB_OSMOL_LMM_3 <- lme4::lmer(data = dab_all,
                    Plasma_Osmol_Rep_mean ~ 
                      Treatment + 
                      Mother_osml +
                      Mass_g + 
                        (1|Mother_ID))

DAB_OSMOL_LMM_3a <- lme4::lmer(data = dab_all,
                    Plasma_Osmol_Rep_mean ~ 
                      Treatment + 
                      Mother_osml +
                      live_neonates +
                        (1|Mother_ID))

summary(DAB_OSMOL_LMM_3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Plasma_Osmol_Rep_mean ~ Treatment + Mother_osml + Mass_g + (1 |  
##     Mother_ID)
##    Data: dab_all
## 
## REML criterion at convergence: 632.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.83638 -0.51605 -0.05075  0.60473  3.11669 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Mother_ID (Intercept)  75.01    8.661  
##  Residual              125.14   11.186  
## Number of obs: 81, groups:  Mother_ID, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 217.2706    49.5555   4.384
## TreatmentW   -8.7918     5.1715  -1.700
## Mother_osml   0.3209     0.1551   2.069
## Mass_g       -1.0816     0.8693  -1.244
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmnW Mthr_s
## TreatmentW  -0.322              
## Mother_osml -0.964  0.324       
## Mass_g      -0.127 -0.200 -0.133
drop1(DAB_OSMOL_LMM_3)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean ~ Treatment + Mother_osml + Mass_g + (1 | 
##     Mother_ID)
##             npar    AIC
## <none>           652.04
## Treatment      1 653.38
## Mother_osml    1 654.78
## Mass_g         1 651.46
anova(DAB_OSMOL_LMM_3, test = "F", type = "2")
## Analysis of Variance Table
##             npar Sum Sq Mean Sq F value
## Treatment      1 916.99  916.99  7.3279
## Mother_osml    1 461.74  461.74  3.6898
## Mass_g         1 193.73  193.73  1.5481
anova(DAB_OSMOL_LMM_3a, test = "F", type = "2")
## Analysis of Variance Table
##               npar Sum Sq Mean Sq F value
## Treatment        1 995.25  995.25  7.7569
## Mother_osml      1 502.92  502.92  3.9197
## live_neonates    1 109.49  109.49  0.8534
DAB_OSMOL_LMM_4 <- lme4::lmer(data = dab_all,
                    Plasma_Osmol_Rep_mean ~ 
                      Treatment + 
                      Mother_osml +
                        (1|Mother_ID))

drop1(DAB_OSMOL_LMM_4)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean ~ Treatment + Mother_osml + (1 | Mother_ID)
##             npar    AIC
## <none>           651.46
## Treatment      1 653.97
## Mother_osml    1 653.70
car::Anova(DAB_OSMOL_LMM_4, test = "F", type = "2")
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: Plasma_Osmol_Rep_mean
##                  F Df Df.res  Pr(>F)  
## Treatment   4.2611  1 14.764 0.05703 .
## Mother_osml 3.9788  1 15.258 0.06426 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DAB_OSMOL_LMM_5 <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      Treatment +
                        (1|Mother_ID))
DAB_OSMOL_LMM_null <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 1 + (1|Mother_ID))

Model Selection

osml_models_all <- list(DAB_OSMOL_LMM_1c, DAB_OSMOL_LMM_2,
                        DAB_OSMOL_LMM_3, DAB_OSMOL_LMM_3a, 
                        DAB_OSMOL_LMM_4, DAB_OSMOL_LMM_5,
                        DAB_OSMOL_LMM_null)

#specify model names
osml_mod_names_all <- c('model1', 'model2',  
                    'model3', 'model3a', 'model4', 'model5', 
                    'modelnull')

#calculate AIC of each model
osml_AICc_all <- data.frame(aictab(cand.set = osml_models_all, 
                                 modnames = osml_mod_names_all))
osml_AICc_all
##    Modnames K     AICc Delta_AICc    ModelLik      AICcWt    Res.LL    Cum.Wt
## 2    model2 7 643.6723  0.0000000 1.000000000 0.309640917 -314.0690 0.3096409
## 1    model1 8 644.3037  0.6313507 0.729296184 0.225819939 -313.1518 0.5354609
## 3    model3 6 645.4118  1.7395217 0.419051753 0.129755569 -316.1384 0.6652164
## 4   model3a 6 645.5000  1.8276648 0.400984542 0.124161221 -316.1824 0.7893776
## 6    model5 4 645.6023  1.9299350 0.380995582 0.117971821 -318.5380 0.9073495
## 5    model4 5 646.1160  2.4436508 0.294691750 0.091248624 -317.6580 0.9985981
## 7 modelnull 3 654.4675 10.7951590 0.004527526 0.001401907 -324.0779 1.0000000

They are literally all equivalent… lol

Tmt only

DAB_OSMOL_tmt_LMM_1 <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      Treatment * 
                        (Mother_Days_in_Treatment + live_neonates) + 
                        Mass_g +
                        (1|Mother_ID))

car::vif(DAB_OSMOL_tmt_LMM_1)                     
##                          Treatment           Mother_Days_in_Treatment 
##                          20.882880                           1.748671 
##                      live_neonates                             Mass_g 
##                           2.240056                           1.394589 
## Treatment:Mother_Days_in_Treatment            Treatment:live_neonates 
##                           6.147350                          22.128174
drop1(DAB_OSMOL_tmt_LMM_1)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean ~ Treatment * (Mother_Days_in_Treatment + 
##     live_neonates) + Mass_g + (1 | Mother_ID)
##                                    npar    AIC
## <none>                                  651.14
## Mass_g                                1 651.31
## Treatment:Mother_Days_in_Treatment    1 649.49
## Treatment:live_neonates               1 653.30
DAB_OSMOL_tmt_LMM_1a <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      Treatment * 
                        Mother_Days_in_Treatment + 
                        live_neonates + 
                        Mass_g +
                        (1|Mother_ID))

car::vif(DAB_OSMOL_tmt_LMM_1a)                     
##                          Treatment           Mother_Days_in_Treatment 
##                           4.638489                           1.744984 
##                      live_neonates                             Mass_g 
##                           1.293482                           1.220707 
## Treatment:Mother_Days_in_Treatment 
##                           5.676306
drop1(DAB_OSMOL_tmt_LMM_1a)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean ~ Treatment * Mother_Days_in_Treatment + 
##     live_neonates + Mass_g + (1 | Mother_ID)
##                                    npar    AIC
## <none>                                  653.30
## live_neonates                         1 652.55
## Mass_g                                1 651.96
## Treatment:Mother_Days_in_Treatment    1 652.50
DAB_OSMOL_tmt_LMM_2a <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      Treatment * 
                        Mother_Days_in_Treatment + 
                        live_neonates + 
                        (1|Mother_ID))

car::vif(DAB_OSMOL_tmt_LMM_2a)                     
##                          Treatment           Mother_Days_in_Treatment 
##                           4.502250                           1.688108 
##                      live_neonates Treatment:Mother_Days_in_Treatment 
##                           1.118255                           5.291212
drop1(DAB_OSMOL_tmt_LMM_2a)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean ~ Treatment * Mother_Days_in_Treatment + 
##     live_neonates + (1 | Mother_ID)
##                                    npar    AIC
## <none>                                  651.96
## live_neonates                         1 650.72
## Treatment:Mother_Days_in_Treatment    1 651.93
DAB_OSMOL_tmt_LMM_2b <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      Treatment + 
                        Mother_Days_in_Treatment + 
                        live_neonates + 
                        Mass_g +
                        (1|Mother_ID))

car::vif(DAB_OSMOL_tmt_LMM_2b)                     
##                Treatment Mother_Days_in_Treatment            live_neonates 
##                 1.035046                 1.019044                 1.135896 
##                   Mass_g 
##                 1.139650
drop1(DAB_OSMOL_tmt_LMM_2b)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean ~ Treatment + Mother_Days_in_Treatment + 
##     live_neonates + Mass_g + (1 | Mother_ID)
##                          npar    AIC
## <none>                        652.50
## Treatment                   1 658.21
## Mother_Days_in_Treatment    1 655.31
## live_neonates               1 652.93
## Mass_g                      1 651.93
DAB_OSMOL_tmt_LMM_3a <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      Treatment * 
                        Mother_Days_in_Treatment + 
                        (1|Mother_ID))

drop1(DAB_OSMOL_tmt_LMM_3a)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean ~ Treatment * Mother_Days_in_Treatment + 
##     (1 | Mother_ID)
##                                    npar    AIC
## <none>                                  650.72
## Treatment:Mother_Days_in_Treatment    1 651.49
DAB_OSMOL_tmt_LMM_3b <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      Treatment + 
                        Mother_Days_in_Treatment + 
                        live_neonates + 
                        (1|Mother_ID))

drop1(DAB_OSMOL_tmt_LMM_3b)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean ~ Treatment + Mother_Days_in_Treatment + 
##     live_neonates + (1 | Mother_ID)
##                          npar    AIC
## <none>                        651.93
## Treatment                   1 658.90
## Mother_Days_in_Treatment    1 654.97
## live_neonates               1 651.49
DAB_OSMOL_tmt_LMM_4 <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      Treatment +
                        Mother_Days_in_Treatment + 
                        (1|Mother_ID))

drop1(DAB_OSMOL_tmt_LMM_4)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean ~ Treatment + Mother_Days_in_Treatment + 
##     (1 | Mother_ID)
##                          npar    AIC
## <none>                        651.49
## Treatment                   1 657.43
## Mother_Days_in_Treatment    1 653.70
anova(DAB_OSMOL_tmt_LMM_4, test = "F", type = "2")
## Analysis of Variance Table
##                          npar  Sum Sq Mean Sq F value
## Treatment                   1 1010.05 1010.05  7.8605
## Mother_Days_in_Treatment    1  507.32  507.32  3.9481
DAB_OSMOL_tmt_LMM_5 <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      Treatment +
                        (1|Mother_ID))

DAB_OSMOL_tmt_LMM_null <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      1 +
                        (1|Mother_ID))

Model Selection

osml_tmt_models_all <- list(DAB_OSMOL_tmt_LMM_1a, 
                            DAB_OSMOL_tmt_LMM_2a, DAB_OSMOL_tmt_LMM_2b,
                        DAB_OSMOL_tmt_LMM_3a, DAB_OSMOL_tmt_LMM_3b, 
                        DAB_OSMOL_tmt_LMM_4, DAB_OSMOL_tmt_LMM_5,
                        DAB_OSMOL_tmt_LMM_null)

#specify model names
osml_tmt_mod_names_all <- c('model1', 'model2a', 'model2b',    
                    'model3a', 'model3b', 'model4', 'model5', 
                    'modelnull')

#calculate AIC of each model
data.frame(aictab(cand.set = osml_tmt_models_all, 
                  modnames = osml_tmt_mod_names_all))
##    Modnames K     AICc Delta_AICc    ModelLik       AICcWt    Res.LL    Cum.Wt
## 1    model1 8 641.8093  0.0000000 1.000000000 0.2108152701 -311.9047 0.2108153
## 2   model2a 7 641.9359  0.1266004 0.938661650 0.1978842093 -313.2009 0.4086995
## 3   model2b 7 642.0460  0.2366023 0.888428466 0.1872942871 -313.2559 0.5959938
## 4   model3a 6 642.2566  0.4472416 0.799618302 0.1685717483 -314.5607 0.7645655
## 5   model3b 6 642.8059  0.9965279 0.607584552 0.1280881015 -314.8354 0.8926536
## 6    model4 5 643.8676  2.0582860 0.357313048 0.0753270467 -316.5338 0.9679807
## 7    model5 4 645.6023  3.7929112 0.150099690 0.0316433066 -318.5380 0.9996240
## 8 modelnull 3 654.4675 12.6581353 0.001783696 0.0003760304 -324.0779 1.0000000

Mom Osml only

DAB_OSMOL_mom_LMM_1 <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      Mother_osml +
                        Mother_Days_in_Treatment + live_neonates + 
                        Mass_g +
                        (1|Mother_ID))

car::vif(DAB_OSMOL_mom_LMM_1)                     
##              Mother_osml Mother_Days_in_Treatment            live_neonates 
##                 1.281577                 1.286981                 1.124035 
##                   Mass_g 
##                 1.122335
drop1(DAB_OSMOL_mom_LMM_1)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean ~ Mother_osml + Mother_Days_in_Treatment + 
##     live_neonates + Mass_g + (1 | Mother_ID)
##                          npar    AIC
## <none>                        655.10
## Mother_osml                 1 658.21
## Mother_Days_in_Treatment    1 653.56
## live_neonates               1 655.09
## Mass_g                      1 656.84
DAB_OSMOL_mom_LMM_2 <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      Mother_osml +
                      live_neonates + 
                      Mass_g +
                        (1|Mother_ID))

drop1(DAB_OSMOL_mom_LMM_2)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean ~ Mother_osml + live_neonates + Mass_g + 
##     (1 | Mother_ID)
##               npar    AIC
## <none>             653.56
## Mother_osml      1 659.09
## live_neonates    1 653.38
## Mass_g           1 655.46
anova(DAB_OSMOL_mom_LMM_2, test = "F", type = "2")
## Analysis of Variance Table
##               npar Sum Sq Mean Sq F value
## Mother_osml      1 782.58  782.58  6.3179
## live_neonates    1  56.14   56.14  0.4533
## Mass_g           1 464.00  464.00  3.7460
DAB_OSMOL_mom_LMM_3 <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      Mother_osml +
                      Mass_g +
                        (1|Mother_ID))

drop1(DAB_OSMOL_mom_LMM_3)
## Single term deletions
## 
## Model:
## Plasma_Osmol_Rep_mean ~ Mother_osml + Mass_g + (1 | Mother_ID)
##             npar    AIC
## <none>           653.38
## Mother_osml    1 658.03
## Mass_g         1 653.97
DAB_OSMOL_mom_LMM_4 <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      Mother_osml +
                        (1|Mother_ID))
anova(DAB_OSMOL_mom_LMM_4, test = "F", type = "2")
## Analysis of Variance Table
##             npar Sum Sq Mean Sq F value
## Mother_osml    1 824.95  824.95  6.4084
DAB_OSMOL_mom_LMM_null <- lme4::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      1 +
                        (1|Mother_ID))

Model Selection

osml_mom_models_all <- list(DAB_OSMOL_mom_LMM_1, 
                            DAB_OSMOL_mom_LMM_2, 
                            DAB_OSMOL_mom_LMM_3, 
                            DAB_OSMOL_mom_LMM_4, 
                            DAB_OSMOL_mom_LMM_null)

#specify model names
osml_mom_mod_names_all <- c('model1', 'model2', 
                            'model3', 'model4', 
                            'modelnull')

#calculate AIC of each model
data.frame(aictab(cand.set = osml_mom_models_all, 
                  modnames = osml_mom_mod_names_all))
##    Modnames K     AICc Delta_AICc   ModelLik     AICcWt    Res.LL    Cum.Wt
## 2    model2 6 649.3303   0.000000 1.00000000 0.45058843 -318.0976 0.4505884
## 1    model1 7 650.6159   1.285680 0.52579701 0.23691805 -317.5408 0.6875065
## 3    model3 5 650.9535   1.623237 0.44413854 0.20012369 -320.0767 0.8876302
## 4    model4 4 652.8422   3.511903 0.17274283 0.07783592 -322.1579 0.9654661
## 5 modelnull 3 654.4675   5.137226 0.07664179 0.03453390 -324.0779 1.0000000

Final Modelsss…

Okay, so if we use either tmt or mother osml, then the top models are one with Treatment * Mother_Days_in_Treatment and one with only Mother_osml!

DAB_OSMOL_LMM_BEST_tmt <- lmerTest::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      Treatment * Mother_Days_in_Treatment + 
                        (1|Mother_ID))
DAB_OSMOL_LMM_BEST_mom <- lmerTest::lmer(data = dab_all,
                      Plasma_Osmol_Rep_mean ~ 
                      Mother_osml +
                        (1|Mother_ID))

anova(DAB_OSMOL_LMM_BEST_tmt, test = "F", type = "2")
## Type II Analysis of Variance Table with Satterthwaite's method
##                                    Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)
## Treatment                            2.49    2.49     1 13.444  0.0194 0.89119
## Mother_Days_in_Treatment           534.04  534.04     1 14.319  4.1664 0.06011
## Treatment:Mother_Days_in_Treatment 301.39  301.39     1 14.248  2.3514 0.14707
##                                     
## Treatment                           
## Mother_Days_in_Treatment           .
## Treatment:Mother_Days_in_Treatment  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(DAB_OSMOL_LMM_BEST_mom, test = "F", type = "2")
## Type II Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)  
## Mother_osml 824.95  824.95     1 16.028  6.4084 0.0222 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans

summary(DAB_OSMOL_LMM_BEST_tmt)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Plasma_Osmol_Rep_mean ~ Treatment * Mother_Days_in_Treatment +  
##     (1 | Mother_ID)
##    Data: dab_all
## 
## REML criterion at convergence: 629.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9034 -0.5272 -0.1034  0.6063  3.1029 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Mother_ID (Intercept)  59.6     7.72   
##  Residual              128.2    11.32   
## Number of obs: 81, groups:  Mother_ID, 18
## 
## Fixed effects:
##                                     Estimate Std. Error       df t value
## (Intercept)                         291.0394     5.9030  13.2003  49.303
## TreatmentW                           -1.2598     9.0369  13.4438  -0.139
## Mother_Days_in_Treatment              1.4761     0.5751  14.4829   2.567
## TreatmentW:Mother_Days_in_Treatment  -1.3886     0.9056  14.2483  -1.533
##                                     Pr(>|t|)    
## (Intercept)                         2.34e-16 ***
## TreatmentW                            0.8912    
## Mother_Days_in_Treatment              0.0219 *  
## TreatmentW:Mother_Days_in_Treatment   0.1471    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmnW M_D__T
## TreatmentW  -0.653              
## Mthr_Dys__T -0.824  0.538       
## TrtW:M_D__T  0.523 -0.869 -0.635
emmeans(DAB_OSMOL_LMM_BEST_tmt, 
        pairwise ~ Treatment | Mother_Days_in_Treatment)
## $emmeans
## Mother_Days_in_Treatment = 8.33:
##  Treatment emmean   SE   df lower.CL upper.CL
##  C            303 3.35 14.2      296      311
##  W            291 2.98 13.7      284      297
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
## Mother_Days_in_Treatment = 8.33:
##  contrast estimate   SE df t.ratio p.value
##  C - W        12.8 4.48 14   2.864  0.0125
## 
## Degrees-of-freedom method: kenward-roger
osml_tmt_trends <- emtrends(DAB_OSMOL_LMM_BEST_tmt, 
         var = "Mother_Days_in_Treatment", "Treatment")
tidy(osml_tmt_trends)
## # A tibble: 2 × 6
##   Treatment Mother_Days_in_Treatment.trend std.error    df statistic p.value
##   <chr>                              <dbl>     <dbl> <dbl>     <dbl>   <dbl>
## 1 C                                 1.48       0.575  14.4     2.57   0.0220
## 2 W                                 0.0875     0.700  14.0     0.125  0.902
pairs(osml_tmt_trends)
##  contrast estimate    SE   df t.ratio p.value
##  C - W        1.39 0.906 14.2   1.533  0.1472
## 
## Degrees-of-freedom method: kenward-roger
DAB_OSMOL_emmeans_tmt <- emmeans(DAB_OSMOL_LMM_BEST_tmt, 
                      "Treatment", by = "Mother_Days_in_Treatment")
DAB_OSMOL_emmeans_CI <- tidy(confint(DAB_OSMOL_emmeans_tmt))
summary(DAB_OSMOL_LMM_BEST_mom)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Plasma_Osmol_Rep_mean ~ Mother_osml + (1 | Mother_ID)
##    Data: dab_all
## 
## REML criterion at convergence: 644.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6944 -0.4909 -0.1177  0.5547  3.2989 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Mother_ID (Intercept)  85.15    9.228  
##  Residual              128.73   11.346  
## Number of obs: 81, groups:  Mother_ID, 18
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)   
## (Intercept) 174.3195    48.2604  15.9829   3.612  0.00234 **
## Mother_osml   0.3895     0.1538  16.0276   2.531  0.02220 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Mother_osml -0.999
r.squaredGLMM(DAB_OSMOL_LMM_BEST_mom)
##           R2m       R2c
## [1,] 0.155193 0.4915269
r.squaredGLMM(DAB_OSMOL_LMM_BEST_tmt)
##            R2m       R2c
## [1,] 0.2675243 0.5000087

Model Assumptions

plot(DAB_OSMOL_LMM_BEST_tmt)

shapiro.test(residuals(DAB_OSMOL_LMM_BEST_tmt))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(DAB_OSMOL_LMM_BEST_tmt)
## W = 0.97279, p-value = 0.08199
plot(DAB_OSMOL_LMM_BEST_mom)

shapiro.test(residuals(DAB_OSMOL_LMM_BEST_mom))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(DAB_OSMOL_LMM_BEST_mom)
## W = 0.97018, p-value = 0.05555

looks great

Partial Regression

actually, I do not need to get partreg values because these relationships are simple/straightforward to plot

plot(effect("Mother_osml", DAB_OSMOL_LMM_BEST_mom))

plot(effect("Treatment:Mother_Days_in_Treatment", 
            DAB_OSMOL_LMM_BEST_tmt), multiline = T)

PUB FIG

osml ~ mom osml

ggplot(
  dab_clutch_mom,
  aes(
    x = Mother_osml,
    y = Plasma_Osmol_Rep_mean_mean
  )
) +
  geom_abline(
    slope = 1, intercept = 0,
    linetype = "dashed", color = "darkgrey"
  ) +
  geom_pointrange(
    aes(ymin = Plasma_Osmol_Rep_mean_mean-Plasma_Osmol_Rep_mean_sd,
        ymax = Plasma_Osmol_Rep_mean_mean+Plasma_Osmol_Rep_mean_sd,
        color = Treatment, shape = Treatment),
    alpha = 0.6, size = 0.2, linewidth = 0.3, show.legend = FALSE
  ) +
  geom_point(aes(color = Treatment, shape = Treatment)) +
  geom_smooth(
    se = T,
    formula = 'y~x',
    method = "lm",
    linewidth = 1,
    color = "black",
    show.legend = FALSE
  ) +
  
  # prettys 
  scale_color_manual(
    values = c("darkgoldenrod2", "royalblue4"), 
    labels = c("C" = "Control", "W" = "Hydration"), 
    name = NULL
  ) +
  scale_shape_manual(
    values = c(15, 17), 
    labels = c("C" = "Control", "W" = "Hydration"),
    name = NULL
  ) +
  theme(
    legend.position = "inside",
    legend.position.inside = c(0.75, 0.15),
    legend.background = element_blank()
  ) +
  
  scale_x_continuous(
    limits = c(260, 360),
    breaks = c(260, 280, 300, 320, 340, 360),
    name = expression('Mother Osmolality (mmol '*kg^-1*')'),
  ) +
  scale_y_continuous(
    limits = c(260, 360),
    breaks = c(260, 280, 300, 320, 340, 360),
    name = expression('Neonate Osmolality (mmol '*kg^-1*')'),
  ) -> neo_mom_osml
neo_mom_osml

export it:

ggsave(
  filename = "Fig5_neonate_mother_osml.pdf",
  plot = neo_mom_osml,
  path = "./results_figures",
  device = "pdf",
  dpi = 600,
  units = "mm",
  width = 150,
  height = 100
)
ggsave(
  filename = "Fig5_neonate_mother_osml.tiff",
  plot = neo_mom_osml,
  path = "./results_figures",
  device = "tiff",
  dpi = 600,
  units = "mm",
  width = 150,
  height = 100
)

osml ~ tmt * days

ggplot(
  dab_clutch_mom,
  aes(
    x = Mother_Days_in_Treatment,
    y = Plasma_Osmol_Rep_mean_mean,
    color = Treatment, 
    shape = Treatment
  )
) +
  # geom_abline(
  #   slope = 1, intercept = 0, 
  #   linetype = "dashed", color = "darkgrey"
  # ) +
  geom_pointrange(
    aes(ymin = Plasma_Osmol_Rep_mean_mean-Plasma_Osmol_Rep_mean_sd,
        ymax = Plasma_Osmol_Rep_mean_mean+Plasma_Osmol_Rep_mean_sd),
    alpha = 0.6, size = 0.2, linewidth = 0.3, show.legend = FALSE
  ) +
  geom_point() +
  geom_smooth(
    se = F,
    formula = 'y~x',
    method = "lm",
    linewidth = 1,
    show.legend = FALSE
  ) +
  
  # prettys 
  scale_color_manual(
    values = c("darkgoldenrod2", "royalblue4"), 
    labels = c("C" = "Control", "W" = "Hydration"), 
    name = NULL
  ) +
  scale_shape_manual(
    values = c(15, 17), 
    labels = c("C" = "Control", "W" = "Hydration"),
    name = NULL
  ) +
  theme(
    legend.position = "inside",
    legend.position.inside = c(0.15, 0.85),
    legend.background = element_blank()
  ) +
  
  scale_x_continuous(
    limits = c(2, 15),
    breaks = c(2, 6, 10, 14),
    name = "Mother Days in Treatment Before Birthing"
  ) +
  scale_y_continuous(
    limits = c(269, 360),
    breaks = c(280, 300, 320, 340, 360),
    name = expression('Neonate Osmolality (mmol '*kg^-1*')'),
  ) -> neo_tmt_osml
neo_tmt_osml

PUB FIG 3 DAB

Osmolality

ggplot() +
  geom_pointrange(
    data = dab_all_clutch_avg,
    aes(
      x = Treatment,
      y = Plasma_Osmol_Rep_mean_mean,
      ymin = Plasma_Osmol_Rep_mean_min,
      ymax = Plasma_Osmol_Rep_mean_max,
      group = Mother_ID,
      color = Treatment,
      shape = Treatment
    ),
    position = position_jitter(width = 0.3),
    size = 0.2, linewidth = 0.3,
    alpha = 0.8,
    show.legend = FALSE
  ) +
  geom_errorbar(
    data = DAB_OSMOL_emmeans_CI,
    aes(x = Treatment, ymin = conf.low, ymax = conf.high,
        color = Treatment),
    width = 0.3
  ) +
  geom_point(
    data = DAB_OSMOL_emmeans_CI,
    aes(x = Treatment, y = estimate, 
        color = Treatment, shape = Treatment),
    size = 2.5
  ) +

  # prettys 
  scale_color_manual(
    values = c("darkgoldenrod2", "royalblue4"), 
    labels = c("C" = "Control", "W" = "Hydration"), 
    name = NULL
  ) +
  scale_shape_manual(
    values = c(15, 17), 
    labels = c("C" = "Control", "W" = "Hydration"),
    name = NULL
  ) +
  labs(
    x = NULL, 
    y = expression('Osmolality (mmol '*kg^-1*')'),
    color = "Treatment"
  ) +
  scale_y_continuous(
    limits = c(260, 360), 
    breaks = seq(260, 360, by = 20)
  ) +
  scale_x_discrete(labels = c("C" = "Control", "W" = "Hydration")) +
  
  theme(legend.position = "none") +
  
  # statistical labels
  annotate(geom = "text", x = 1.5, y = 340, label = "*", size = 6) -> dab_osml
  
dab_osml

CEWL

ggplot() +
  geom_pointrange(
    data = dab_all_clutch_avg,
    aes(
      x = Treatment,
      y = CEWL_g_m2h_mean,
      ymin = CEWL_g_m2h_min,
      ymax = CEWL_g_m2h_max,
      group = Mother_ID,
      color = Treatment,
      shape = Treatment
    ),
    position = position_jitter(width = 0.3),
    size = 0.2, linewidth = 0.3,
    alpha = 0.8,
    show.legend = FALSE
  ) +
  geom_errorbar(
    data = DAB_CEWL_emmeans_CI,
    aes(x = Treatment, ymin = conf.low, ymax = conf.high,
        color = Treatment),
    width = 0.3
  ) +
  geom_point(
    data = DAB_CEWL_emmeans_CI,
    aes(x = Treatment, y = estimate, 
        color = Treatment, shape = Treatment),
    size = 2.5
  ) +

  # prettys 
  scale_color_manual(
    values = c("darkgoldenrod2", "royalblue4"), 
    labels = c("C" = "Control", "W" = "Hydration"), 
    name = NULL
  ) +
  scale_shape_manual(
    values = c(15, 17), 
    labels = c("C" = "Control", "W" = "Hydration"),
    name = NULL
  ) +
  labs(
    x = NULL, 
    y = expression('CEWL (g '*m^-2*' '*h^-1*')'), 
    color = "Treatment"
  ) +
  scale_y_continuous(
    limits = c(4, 25), 
    breaks = seq(0, 30, by = 5)
  ) +
  scale_x_discrete(labels = c("C" = "Control", "W" = "Hydration")) +
  
  theme(legend.position = "none") +
  
  # statistical labels
  annotate(geom = "text", x = 1.5, y = 22, label = "NS", size = 2) -> dab_cewl
  
dab_cewl

Arrange

set.seed(135) # to get a nice jitter: 10, 35, 44, 359, 135

ggarrange(
  dab_osml, dab_cewl,
  nrow = 2,
  ncol = 1,
  labels = c("A", "B"),
  align = "hv"
) -> dab_diffs

ggsave(
  filename = "Fig3_neonate_dab_diffs.pdf",
  plot = dab_diffs,
  path = "./results_figures",
  device = "pdf",
  dpi = 600,
  units = "mm",
  width = 80,
  height = 120
)

ggsave(
  filename = "Fig3_neonate_dab_diffs.tiff",
  plot = dab_diffs,
  path = "./results_figures",
  device = "tiff",
  dpi = 600,
  units = "mm",
  width = 80,
  height = 120
)

TMT EFFECTS

Calculate Deltas

Only 3 individuals were for sure measured at the before and after time points, so the change from before to after cannot be calculated for individual neonates.

deltas_means <- dab_ps_compare_means %>%
  arrange(Mother_ID, timept) %>% 
  group_by(Mother_ID) %>%
  mutate(
    delta_mass = Mass_g_mean - lag(Mass_g_mean),
    delta_osml = Plasma_Osmol_Rep_mean_mean - lag(Plasma_Osmol_Rep_mean_mean),
    delta_CEWL = CEWL_g_m2h_mean - lag(CEWL_g_m2h_mean)
  ) 

deltas_befaft <- deltas_means %>% 
  ungroup() %>% 
  mutate(
    delta_mass = case_when(is.na(delta_mass) ~ 0, TRUE ~ delta_mass),
    delta_osml = case_when(is.na(delta_osml) ~ 0, TRUE ~ delta_osml),
    delta_CEWL = case_when(is.na(delta_CEWL) ~ 0, TRUE ~ delta_CEWL)
  )

deltas <- deltas_means %>%
  filter(complete.cases(delta_mass))

summary(deltas_befaft)
##    Mother_ID Treatment timept        n         Tail_Length_cm_mean
##  124    :2   C:8       DAB:8   Min.   :3.000   Min.   :1.560      
##  125    :2   W:8       PS :8   1st Qu.:4.750   1st Qu.:1.775      
##  128    :2                     Median :5.000   Median :1.880      
##  129    :2                     Mean   :4.625   Mean   :1.865      
##  131    :2                     3rd Qu.:5.000   3rd Qu.:1.962      
##  133    :2                     Max.   :5.000   Max.   :2.160      
##  (Other):4                                                        
##   SVL_cm_mean     Mass_g_mean    Tb_CEWL_c_mean  Plasma_Osmol_Rep_mean_mean
##  Min.   :21.80   Min.   :11.66   Min.   :27.42   Min.   :280.5             
##  1st Qu.:22.31   1st Qu.:13.07   1st Qu.:27.96   1st Qu.:286.2             
##  Median :23.35   Median :14.14   Median :28.60   Median :296.6             
##  Mean   :23.17   Mean   :14.40   Mean   :28.58   Mean   :304.2             
##  3rd Qu.:23.72   3rd Qu.:15.07   3rd Qu.:29.23   3rd Qu.:312.5             
##  Max.   :25.87   Max.   :18.46   Max.   :29.92   Max.   :363.4             
##                                                                            
##  CEWL_g_m2h_mean  msmt_temp_C_mean msmt_RH_percent_mean  VPD_kPa_mean  
##  Min.   : 7.296   Min.   :23.71    Min.   :36.80        Min.   :1.595  
##  1st Qu.:10.795   1st Qu.:25.40    1st Qu.:39.02        1st Qu.:1.923  
##  Median :15.046   Median :25.75    Median :40.10        Median :1.962  
##  Mean   :13.959   Mean   :25.65    Mean   :40.39        Mean   :1.965  
##  3rd Qu.:16.193   3rd Qu.:26.01    3rd Qu.:41.24        3rd Qu.:2.060  
##  Max.   :17.244   Max.   :26.88    Max.   :45.57        Max.   :2.160  
##                                                                        
##  Tail_Length_cm_sd   SVL_cm_sd        Mass_g_sd       Tb_CEWL_c_sd   
##  Min.   :0.1517    Min.   :0.4435   Min.   :0.2380   Min.   :0.3000  
##  1st Qu.:0.1930    1st Qu.:0.6132   1st Qu.:0.4584   1st Qu.:0.4959  
##  Median :0.2334    Median :0.7606   Median :0.6737   Median :0.9314  
##  Mean   :0.2407    Mean   :0.7866   Mean   :0.8755   Mean   :0.9654  
##  3rd Qu.:0.2842    3rd Qu.:0.9680   3rd Qu.:1.2115   3rd Qu.:1.3562  
##  Max.   :0.3808    Max.   :1.0970   Max.   :2.0768   Max.   :1.8050  
##                                                                      
##  Plasma_Osmol_Rep_mean_sd CEWL_g_m2h_sd   msmt_temp_C_sd   msmt_RH_percent_sd
##  Min.   : 2.739           Min.   :1.128   Min.   :0.1311   Min.   :0.1248    
##  1st Qu.: 6.355           1st Qu.:1.765   1st Qu.:0.1642   1st Qu.:0.3391    
##  Median : 8.009           Median :2.213   Median :0.1990   Median :0.6041    
##  Mean   :12.556           Mean   :2.577   Mean   :0.2269   Mean   :0.8738    
##  3rd Qu.:14.390           3rd Qu.:3.172   3rd Qu.:0.2876   3rd Qu.:1.2300    
##  Max.   :43.634           Max.   :5.040   Max.   :0.3929   Max.   :2.6074    
##                                                                              
##    VPD_kPa_sd       Tail_Length_cm_min   SVL_cm_min      Mass_g_min   
##  Min.   :0.009671   Min.   :1.200      Min.   :20.50   Min.   :11.10  
##  1st Qu.:0.028434   1st Qu.:1.500      1st Qu.:21.38   1st Qu.:11.90  
##  Median :0.037982   Median :1.500      Median :22.45   Median :13.05  
##  Mean   :0.044378   Mean   :1.562      Mean   :22.18   Mean   :13.23  
##  3rd Qu.:0.062645   3rd Qu.:1.700      3rd Qu.:23.00   3rd Qu.:14.28  
##  Max.   :0.110722   Max.   :1.900      Max.   :24.60   Max.   :16.60  
##                                                                       
##  Tb_CEWL_c_min   Plasma_Osmol_Rep_mean_min CEWL_g_m2h_min   msmt_temp_C_min
##  Min.   :25.80   Min.   :242.7             Min.   : 4.745   Min.   :23.56  
##  1st Qu.:26.98   1st Qu.:279.5             1st Qu.: 9.571   1st Qu.:25.06  
##  Median :27.40   Median :286.2             Median :11.494   Median :25.52  
##  Mean   :27.50   Mean   :290.2             Mean   :10.959   Mean   :25.40  
##  3rd Qu.:27.93   3rd Qu.:296.2             3rd Qu.:12.982   3rd Qu.:25.76  
##  Max.   :29.50   Max.   :332.5             Max.   :14.812   Max.   :26.70  
##                                                                            
##  msmt_RH_percent_min  VPD_kPa_min    Tail_Length_cm_max   SVL_cm_max   
##  Min.   :35.15       Min.   :1.566   Min.   :2.000      Min.   :22.40  
##  1st Qu.:37.66       1st Qu.:1.858   1st Qu.:2.000      1st Qu.:23.45  
##  Median :39.23       Median :1.917   Median :2.050      Median :24.00  
##  Mean   :39.44       Mean   :1.911   Mean   :2.119      Mean   :24.11  
##  3rd Qu.:40.96       3rd Qu.:1.972   3rd Qu.:2.200      3rd Qu.:24.62  
##  Max.   :44.34       Max.   :2.144   Max.   :2.500      Max.   :26.50  
##                                                                        
##    Mass_g_max    Tb_CEWL_c_max   Plasma_Osmol_Rep_mean_max CEWL_g_m2h_max  
##  Min.   :12.10   Min.   :28.10   Min.   :283.5             Min.   : 9.705  
##  1st Qu.:13.70   1st Qu.:29.35   1st Qu.:297.5             1st Qu.:13.229  
##  Median :14.70   Median :30.20   Median :306.2             Median :18.466  
##  Mean   :15.25   Mean   :29.83   Mean   :318.5             Mean   :17.158  
##  3rd Qu.:16.38   3rd Qu.:30.43   3rd Qu.:327.8             3rd Qu.:19.771  
##  Max.   :20.00   Max.   :31.30   Max.   :428.0             Max.   :23.808  
##                                                                            
##  msmt_temp_C_max msmt_RH_percent_max  VPD_kPa_max      delta_mass     
##  Min.   :24.02   Min.   :37.77       Min.   :1.641   Min.   :-2.4800  
##  1st Qu.:25.68   1st Qu.:40.20       1st Qu.:1.963   1st Qu.:-1.4400  
##  Median :25.99   Median :41.31       Median :2.001   Median :-0.4100  
##  Mean   :25.93   Mean   :41.46       Mean   :2.014   Mean   :-0.8094  
##  3rd Qu.:26.29   3rd Qu.:42.48       3rd Qu.:2.138   3rd Qu.: 0.0000  
##  Max.   :27.20   Max.   :46.22       Max.   :2.179   Max.   : 0.0000  
##                                                                       
##    delta_osml       delta_CEWL    
##  Min.   :-18.17   Min.   :-1.215  
##  1st Qu.:  0.00   1st Qu.: 0.000  
##  Median :  0.00   Median : 0.000  
##  Mean   : 10.75   Mean   : 2.193  
##  3rd Qu.: 14.84   3rd Qu.: 5.339  
##  Max.   : 72.97   Max.   : 9.816  
## 
summary(deltas)
##    Mother_ID Treatment timept        n         Tail_Length_cm_mean
##  124    :1   C:4       DAB:0   Min.   :3.000   Min.   :1.680      
##  125    :1   W:4       PS :8   1st Qu.:3.750   1st Qu.:1.775      
##  128    :1                     Median :5.000   Median :1.890      
##  129    :1                     Mean   :4.375   Mean   :1.890      
##  131    :1                     3rd Qu.:5.000   3rd Qu.:1.983      
##  133    :1                     Max.   :5.000   Max.   :2.160      
##  (Other):2                                                        
##   SVL_cm_mean     Mass_g_mean    Tb_CEWL_c_mean  Plasma_Osmol_Rep_mean_mean
##  Min.   :23.33   Min.   :11.66   Min.   :27.47   Min.   :282.2             
##  1st Qu.:23.48   1st Qu.:12.52   1st Qu.:27.75   1st Qu.:289.3             
##  Median :23.64   Median :13.26   Median :28.67   Median :313.0             
##  Mean   :23.94   Mean   :13.59   Mean   :28.54   Mean   :314.9             
##  3rd Qu.:24.07   3rd Qu.:13.80   3rd Qu.:29.28   3rd Qu.:335.1             
##  Max.   :25.87   Max.   :17.53   Max.   :29.54   Max.   :363.4             
##                                                                            
##  CEWL_g_m2h_mean msmt_temp_C_mean msmt_RH_percent_mean  VPD_kPa_mean  
##  Min.   :14.57   Min.   :25.07    Min.   :36.80        Min.   :1.924  
##  1st Qu.:15.57   1st Qu.:25.49    1st Qu.:37.87        1st Qu.:1.953  
##  Median :16.34   Median :25.62    Median :39.35        Median :1.967  
##  Mean   :16.15   Mean   :25.67    Mean   :39.21        Mean   :2.003  
##  3rd Qu.:17.07   3rd Qu.:26.00    3rd Qu.:40.12        3rd Qu.:2.060  
##  Max.   :17.24   Max.   :26.02    Max.   :41.79        Max.   :2.124  
##                                                                       
##  Tail_Length_cm_sd   SVL_cm_sd        Mass_g_sd       Tb_CEWL_c_sd   
##  Min.   :0.1517    Min.   :0.4435   Min.   :0.2775   Min.   :0.3000  
##  1st Qu.:0.1930    1st Qu.:0.5511   1st Qu.:0.4584   1st Qu.:0.4632  
##  Median :0.2121    Median :0.6728   Median :0.6055   Median :0.7047  
##  Mean   :0.2143    Mean   :0.7485   Mean   :0.8381   Mean   :0.7398  
##  3rd Qu.:0.2255    3rd Qu.:0.9966   3rd Qu.:0.9867   3rd Qu.:0.9812  
##  Max.   :0.2881    Max.   :1.0970   Max.   :1.9655   Max.   :1.2973  
##                                                                      
##  Plasma_Osmol_Rep_mean_sd CEWL_g_m2h_sd   msmt_temp_C_sd   msmt_RH_percent_sd
##  Min.   : 4.174           Min.   :1.595   Min.   :0.1311   Min.   :0.1248    
##  1st Qu.: 7.866           1st Qu.:2.560   1st Qu.:0.1925   1st Qu.:0.3391    
##  Median :11.468           Median :2.967   Median :0.2689   Median :0.5981    
##  Mean   :17.693           Mean   :2.983   Mean   :0.2629   Mean   :0.8022    
##  3rd Qu.:24.097           3rd Qu.:3.565   3rd Qu.:0.3181   3rd Qu.:1.2300    
##  Max.   :43.634           Max.   :4.161   Max.   :0.3929   Max.   :1.7985    
##                                                                              
##    VPD_kPa_sd      Tail_Length_cm_min   SVL_cm_min      Mass_g_min   
##  Min.   :0.02112   Min.   :1.500      Min.   :22.40   Min.   :11.10  
##  1st Qu.:0.03340   1st Qu.:1.500      1st Qu.:22.75   1st Qu.:11.47  
##  Median :0.04012   Median :1.550      Median :23.00   Median :12.15  
##  Mean   :0.04599   Mean   :1.625      Mean   :23.07   Mean   :12.54  
##  3rd Qu.:0.06512   3rd Qu.:1.725      3rd Qu.:23.05   3rd Qu.:13.07  
##  Max.   :0.06985   Max.   :1.900      Max.   :24.60   Max.   :15.30  
##                                                                      
##  Tb_CEWL_c_min   Plasma_Osmol_Rep_mean_min CEWL_g_m2h_min  msmt_temp_C_min
##  Min.   :26.50   Min.   :242.7             Min.   :10.86   Min.   :24.65  
##  1st Qu.:27.23   1st Qu.:283.4             1st Qu.:11.53   1st Qu.:25.27  
##  Median :27.60   Median :296.2             Median :12.28   Median :25.40  
##  Mean   :27.70   Mean   :296.3             Mean   :12.63   Mean   :25.38  
##  3rd Qu.:28.02   3rd Qu.:317.9             3rd Qu.:13.87   3rd Qu.:25.64  
##  Max.   :29.10   Max.   :332.5             Max.   :14.81   Max.   :25.77  
##                                                                           
##  msmt_RH_percent_min  VPD_kPa_min    Tail_Length_cm_max   SVL_cm_max   
##  Min.   :36.20       Min.   :1.859   Min.   :2.000      Min.   :24.00  
##  1st Qu.:37.01       1st Qu.:1.908   1st Qu.:2.000      1st Qu.:24.00  
##  Median :38.16       Median :1.925   Median :2.150      Median :24.25  
##  Mean   :38.48       Mean   :1.945   Mean   :2.125      Mean   :24.77  
##  3rd Qu.:39.74       3rd Qu.:1.969   3rd Qu.:2.200      3rd Qu.:25.55  
##  Max.   :41.52       Max.   :2.063   Max.   :2.300      Max.   :26.50  
##                                                                        
##    Mass_g_max    Tb_CEWL_c_max   Plasma_Osmol_Rep_mean_max CEWL_g_m2h_max 
##  Min.   :12.10   Min.   :28.10   Min.   :296.0             Min.   :18.27  
##  1st Qu.:12.95   1st Qu.:28.27   1st Qu.:300.2             1st Qu.:18.69  
##  Median :14.00   Median :29.85   Median :325.5             Median :19.23  
##  Mean   :14.35   Mean   :29.44   Mean   :334.9             Mean   :19.73  
##  3rd Qu.:14.97   3rd Qu.:30.30   3rd Qu.:352.2             3rd Qu.:20.86  
##  Max.   :19.00   Max.   :30.60   Max.   :428.0             Max.   :21.93  
##                                                                           
##  msmt_temp_C_max msmt_RH_percent_max  VPD_kPa_max      delta_mass    
##  Min.   :25.43   Min.   :37.77       Min.   :1.951   Min.   :-2.480  
##  1st Qu.:25.82   1st Qu.:39.60       1st Qu.:1.996   1st Qu.:-2.235  
##  Median :25.99   Median :40.64       Median :2.009   Median :-1.520  
##  Mean   :25.96   Mean   :40.29       Mean   :2.048   Mean   :-1.619  
##  3rd Qu.:26.14   3rd Qu.:41.18       3rd Qu.:2.138   3rd Qu.:-1.119  
##  Max.   :26.35   Max.   :42.35       Max.   :2.149   Max.   :-0.820  
##                                                                      
##    delta_osml       delta_CEWL    
##  Min.   :-18.17   Min.   :-1.215  
##  1st Qu.:  2.30   1st Qu.: 1.772  
##  Median : 19.85   Median : 5.459  
##  Mean   : 21.50   Mean   : 4.387  
##  3rd Qu.: 35.42   3rd Qu.: 5.942  
##  Max.   : 72.97   Max.   : 9.816  
## 

MASS

mass_lm_time_tmt <- lmerTest::lmer(data = dab_ps_compare, 
            Mass_g ~ timept * Treatment + (1|Mother_ID))
r.squaredGLMM(mass_lm_time_tmt)
##            R2m      R2c
## [1,] 0.1457843 0.817281
mass_means_time <- emmeans(mass_lm_time_tmt, 
                       pairwise ~ timept | Treatment)
mass_means_time$contrasts
## Treatment = C:
##  contrast estimate    SE df t.ratio p.value
##  DAB - PS     2.16 0.319 64   6.760  <.0001
## 
## Treatment = W:
##  contrast estimate    SE df t.ratio p.value
##  DAB - PS     1.08 0.339 64   3.197  0.0022
## 
## Degrees-of-freedom method: kenward-roger
mass_means_tmt <- emmeans(mass_lm_time_tmt, 
                       pairwise ~ Treatment | timept)
mass_means_tmt$contrasts
## timept = DAB:
##  contrast estimate   SE   df t.ratio p.value
##  C - W       0.218 1.38 6.31   0.158  0.8798
## 
## timept = PS:
##  contrast estimate   SE   df t.ratio p.value
##  C - W      -0.853 1.39 6.40  -0.614  0.5606
## 
## Degrees-of-freedom method: kenward-roger
mass_means <- emmeans(mass_lm_time_tmt, "Treatment", by = "timept")
mass_lm_time_tmt_CI <- tidy(confint(mass_means))
#mass_lm_time_tmt_CI <- tidy(mass_means_time$emmeans) # either would work

anova(mass_lm_time_tmt, test = "F", type = "2")
## Type II Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## timept           49.904  49.904     1 64.010 50.4859 1.214e-09 ***
## Treatment         0.043   0.043     1  5.984  0.0440   0.84090    
## timept:Treatment  5.234   5.234     1 64.012  5.2954   0.02465 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

MASS CHANGE

delta_mass_lm <- lm(data = deltas, delta_mass ~ Treatment)

summary(delta_mass_lm)
## 
## Call:
## lm(formula = delta_mass ~ Treatment, data = deltas)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.31500 -0.15812 -0.08292  0.17250  0.48500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.1650     0.1487 -14.560 6.58e-06 ***
## TreatmentW    1.0925     0.2103   5.195  0.00202 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2974 on 6 degrees of freedom
## Multiple R-squared:  0.8181, Adjusted R-squared:  0.7878 
## F-statistic: 26.99 on 1 and 6 DF,  p-value: 0.002024
car::Anova(delta_mass_lm, test = "F", type = "2")
## Anova Table (Type II tests)
## 
## Response: delta_mass
##            Sum Sq Df F value   Pr(>F)   
## Treatment 2.38711  1   26.99 0.002024 **
## Residuals 0.53066  6                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(delta_mass_lm, pairwise ~ Treatment)
## $emmeans
##  Treatment emmean    SE df lower.CL upper.CL
##  C          -2.17 0.149  6    -2.53   -1.801
##  W          -1.07 0.149  6    -1.44   -0.709
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate   SE df t.ratio p.value
##  C - W       -1.09 0.21  6  -5.195  0.0020

OSML

osml_lm_time_tmt <- lmerTest::lmer(data = dab_ps_compare, 
    Plasma_Osmol_Rep_mean ~ timept * Treatment + (1|Mother_ID))
r.squaredGLMM(osml_lm_time_tmt)
##            R2m       R2c
## [1,] 0.6248089 0.6374795
osml_means_time <- emmeans(osml_lm_time_tmt, 
                       pairwise ~ timept | Treatment)
osml_means_time$contrasts
## Treatment = C:
##  contrast estimate   SE   df t.ratio p.value
##  DAB - PS   -44.98 5.30 64.2  -8.484  <.0001
## 
## Treatment = W:
##  contrast estimate   SE   df t.ratio p.value
##  DAB - PS     2.46 5.63 65.0   0.436  0.6644
## 
## Degrees-of-freedom method: kenward-roger
osml_means_tmt <- emmeans(osml_lm_time_tmt, 
                       pairwise ~ Treatment | timept)
osml_means_tmt$contrasts
## timept = DAB:
##  contrast estimate   SE   df t.ratio p.value
##  C - W        5.13 5.73 16.1   0.894  0.3845
## 
## timept = PS:
##  contrast estimate   SE   df t.ratio p.value
##  C - W       52.56 6.06 18.6   8.680  <.0001
## 
## Degrees-of-freedom method: kenward-roger
osml_means <- emmeans(osml_lm_time_tmt, "Treatment", by = "timept")
osml_lm_time_tmt_CI <- tidy(confint(osml_means))

anova(osml_lm_time_tmt, test = "F", type = "2")
## Type II Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## timept            9418.2  9418.2     1 64.945  34.452 1.625e-07 ***
## Treatment        10504.4 10504.4     1  6.297  38.425  0.000675 ***
## timept:Treatment 10315.9 10315.9     1 64.989  37.735 5.492e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

OSML CHANGE

delta_osml_lm <- lm(data = deltas, delta_osml ~ Treatment)

summary(delta_osml_lm)
## 
## Call:
## lm(formula = delta_osml ~ Treatment, data = deltas)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.160 -15.311  -1.113  10.368  26.948 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   46.027      8.664   5.312  0.00181 **
## TreatmentW   -49.062     12.253  -4.004  0.00709 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.33 on 6 degrees of freedom
## Multiple R-squared:  0.7277, Adjusted R-squared:  0.6823 
## F-statistic: 16.03 on 1 and 6 DF,  p-value: 0.007087
car::Anova(delta_osml_lm, test = "F", type = "2")
## Anova Table (Type II tests)
## 
## Response: delta_osml
##           Sum Sq Df F value   Pr(>F)   
## Treatment 4814.2  1  16.032 0.007087 **
## Residuals 1801.7  6                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(delta_osml_lm, pairwise ~ Treatment)
## $emmeans
##  Treatment emmean   SE df lower.CL upper.CL
##  C          46.03 8.66  6     24.8     67.2
##  W          -3.04 8.66  6    -24.2     18.2
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate   SE df t.ratio p.value
##  C - W        49.1 12.3  6   4.004  0.0071

CEWL

NOTE I did test the inclusion of Tb, msmt temp, and msmt VPD in this model specifically, and none had an influence on CEWL compared to the time difference.

cewl_lm_time_tmt <- lmerTest::lmer(data = dab_ps_compare, 
          CEWL_g_m2h ~ timept * Treatment + (1|Mother_ID))
r.squaredGLMM(cewl_lm_time_tmt)
##            R2m       R2c
## [1,] 0.3063975 0.3942069
cewl_means_time <- emmeans(cewl_lm_time_tmt, 
                       pairwise ~ timept | Treatment)
cewl_means_time$contrasts
## Treatment = C:
##  contrast estimate   SE   df t.ratio p.value
##  DAB - PS    -3.67 1.03 63.3  -3.563  0.0007
## 
## Treatment = W:
##  contrast estimate   SE   df t.ratio p.value
##  DAB - PS    -4.64 1.08 63.6  -4.295  0.0001
## 
## Degrees-of-freedom method: kenward-roger
cewl_means_tmt <- emmeans(cewl_lm_time_tmt, 
                       pairwise ~ Treatment | timept)
cewl_means_tmt$contrasts
## timept = DAB:
##  contrast estimate   SE   df t.ratio p.value
##  C - W       1.950 1.34 11.4   1.458  0.1718
## 
## timept = PS:
##  contrast estimate   SE   df t.ratio p.value
##  C - W       0.981 1.38 12.6   0.712  0.4895
## 
## Degrees-of-freedom method: kenward-roger
cewl_means <- emmeans(cewl_lm_time_tmt, "Treatment", by = "timept")
cewl_lm_time_tmt_CI <- tidy(confint(cewl_means))

anova(cewl_lm_time_tmt, test = "F", type = "2")
## Type II Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## timept           309.391 309.391     1 63.603 30.8195 5.913e-07 ***
## Treatment         17.295  17.295     1  6.122  1.7228    0.2364    
## timept:Treatment   4.233   4.233     1 63.615  0.4217    0.5185    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

CEWL CHANGE

delta_cewl_lm <- lm(data = deltas, delta_CEWL ~ Treatment)

summary(delta_cewl_lm)
## 
## Call:
## lm(formula = delta_CEWL ~ Treatment, data = deltas)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2523 -1.9639  0.7369  1.9903  4.7785 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    3.736      1.825   2.047   0.0866 .
## TreatmentW     1.301      2.581   0.504   0.6322  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.651 on 6 degrees of freedom
## Multiple R-squared:  0.04062,    Adjusted R-squared:  -0.1193 
## F-statistic: 0.254 on 1 and 6 DF,  p-value: 0.6322
car::Anova(delta_cewl_lm, test = "F", type = "2")
## Anova Table (Type II tests)
## 
## Response: delta_CEWL
##           Sum Sq Df F value Pr(>F)
## Treatment  3.386  1   0.254 0.6322
## Residuals 79.967  6
emmeans(delta_cewl_lm, pairwise ~ Treatment)
## $emmeans
##  Treatment emmean   SE df lower.CL upper.CL
##  C           3.74 1.83  6   -0.730      8.2
##  W           5.04 1.83  6    0.571      9.5
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate   SE df t.ratio p.value
##  C - W        -1.3 2.58  6  -0.504  0.6322

Tb

tb_lm_time_tmt <- lmerTest::lmer(data = dab_ps_compare, 
          Tb_CEWL_c ~ timept * Treatment + (1|Mother_ID))
r.squaredGLMM(tb_lm_time_tmt)
##             R2m        R2c
## [1,] 0.05646733 0.05646733
tb_means_time <- emmeans(tb_lm_time_tmt, 
                       pairwise ~ timept | Treatment)
tb_means_time$contrasts
## Treatment = C:
##  contrast estimate    SE   df t.ratio p.value
##  DAB - PS   -0.471 0.394 64.3  -1.196  0.2362
## 
## Treatment = W:
##  contrast estimate    SE   df t.ratio p.value
##  DAB - PS    0.486 0.418 65.3   1.161  0.2498
## 
## Degrees-of-freedom method: kenward-roger
tb_means_tmt <- emmeans(tb_lm_time_tmt, 
                       pairwise ~ Treatment | timept)
tb_means_tmt$contrasts
## timept = DAB:
##  contrast estimate    SE   df t.ratio p.value
##  C - W     -0.0921 0.394 19.9  -0.234  0.8175
## 
## timept = PS:
##  contrast estimate    SE   df t.ratio p.value
##  C - W      0.8648 0.419 22.9   2.063  0.0506
## 
## Degrees-of-freedom method: kenward-roger
tb_means <- emmeans(tb_lm_time_tmt, "Treatment", by = "timept")
tb_lm_time_tmt_CI <- tidy(confint(tb_means))

anova(tb_lm_time_tmt, test = "F", type = "2")
## Type II Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## timept           0.0074  0.0074     1    70  0.0049 0.94432  
## Treatment        2.3737  2.3737     1    70  1.5721 0.21407  
## timept:Treatment 4.2049  4.2049     1    70  2.7850 0.09962 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Shedding time

Is there a difference in time between birth and shedding based on tmt group?

date_tmt <- dab_ps_compare %>% 
  group_by(Mother_ID, Treatment) %>% 
  summarise(
    min_date = min(Date_Born_Shed),
    max_date = max(Date_Born_Shed)
  ) %>% 
  mutate(range = as.numeric(max_date - min_date)) %>% 
  arrange(range)

shedtimelm <- lm(data = date_tmt, range ~ Treatment)
summary(shedtimelm)
## 
## Call:
## lm(formula = range ~ Treatment, data = date_tmt)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0000 -0.4375 -0.2500  0.1875  2.0000 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.0000     0.5303   11.31 2.85e-05 ***
## TreatmentW   -0.7500     0.7500   -1.00    0.356    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.061 on 6 degrees of freedom
## Multiple R-squared:  0.1429, Adjusted R-squared:  -2.22e-16 
## F-statistic:     1 on 1 and 6 DF,  p-value: 0.3559
anova(shedtimelm)
## Analysis of Variance Table
## 
## Response: range
##           Df Sum Sq Mean Sq F value Pr(>F)
## Treatment  1  1.125   1.125       1 0.3559
## Residuals  6  6.750   1.125

PUB FIG 6

mass raw

ggplot() +
  geom_pointrange(
    data = dab_ps_compare_means,
    aes(
      x = timept,
      y = Mass_g_mean,
      ymin = Mass_g_mean-Mass_g_sd,
      ymax = Mass_g_mean+Mass_g_sd,
      group = Treatment,
      color = Treatment,
      shape = Treatment
    ),
    position = position_jitterdodge(jitter.width = 0.6, dodge.width = 0.6),
    size = 0.2, linewidth = 0.3,
    alpha = 0.6,
    show.legend = FALSE
  ) +
  geom_line(
    data = mass_lm_time_tmt_CI,
    aes(
      x = timept, y = estimate, 
      group = Treatment, 
      color = Treatment
    ),
    position = position_dodge(width = 0.6)
  ) +
  geom_errorbar(
    data = mass_lm_time_tmt_CI,
    aes(
      x = timept, y = estimate, 
      ymin = conf.low, ymax = conf.high, 
      group = Treatment, 
      color = Treatment
    ),
    width = 0.4,
    position = position_dodge(width = 0.6)
  ) +
  geom_point(
    data = mass_lm_time_tmt_CI,
    aes(
      x = timept, y = estimate, 
      group = Treatment, 
      color = Treatment,
      shape = Treatment
    ),
    size = 2.5,
    position = position_dodge(width = 0.6)
  ) + 
  
  # prettys 
  scale_color_manual(
    values = c("darkgoldenrod2", "royalblue4"), 
    labels = c("C" = "Control", "W" = "Hydration"), 
    name = "Treatment"
  ) +
  scale_shape_manual(
    values = c(15, 17), 
    labels = c("C" = "Control", "W" = "Hydration"),
    name = "Treatment"
  ) +
  scale_y_continuous(
    limits = c(10, 21), 
    breaks = c(10, 15, 20),
    name = "Mass (g)"
  ) + 
  scale_x_discrete(
    labels = c("DAB" = "At Birth", "PS" = "Post-Shed"),
    name = NULL
  ) +
  
  theme(legend.position = "none") +
  
  # statistical labels
  annotate(geom = "text", x = 1, y = 21, label = "NS", size = 2) +
  annotate(geom = "text", x = 2, y = 21, label = "NS", size = 2) +
  annotate(geom = "text", x = 2.5, y = 12.8, label = "*", size = 6, 
           color = "darkgoldenrod2") +
  annotate(geom = "text", x = 2.5, y = 13.7, label = "*", size = 6,
           color = "royalblue4") -> neo_mass # save

neo_mass

osmol raw

ggplot() +
  geom_pointrange(
    data = dab_ps_compare_means,
    aes(
      x = timept,
      y = Plasma_Osmol_Rep_mean_mean,
      ymin = Plasma_Osmol_Rep_mean_mean-Plasma_Osmol_Rep_mean_sd,
      ymax = Plasma_Osmol_Rep_mean_mean+Plasma_Osmol_Rep_mean_sd,
      group = Treatment,
      color = Treatment,
      shape = Treatment
    ),
    position = position_jitterdodge(jitter.width = 0.6, dodge.width = 0.6),
    size = 0.2, linewidth = 0.3,
    alpha = 0.6,
    show.legend = FALSE
  ) +
  geom_line(
    data = osml_lm_time_tmt_CI,
    aes(
      x = timept, y = estimate, 
      group = Treatment, 
      color = Treatment
    ),
    position = position_dodge(width = 0.6)
  ) +
  geom_errorbar(
    data = osml_lm_time_tmt_CI,
    aes(
      x = timept, y = estimate, 
      ymin = conf.low, ymax = conf.high, 
      group = Treatment, 
      color = Treatment
    ),
    width = 0.4,
    position = position_dodge(width = 0.6)
  ) +
  geom_point(
    data = osml_lm_time_tmt_CI,
    aes(
      x = timept, y = estimate, 
      group = Treatment, 
      color = Treatment,
      shape = Treatment
    ),
    size = 2.5,
    position = position_dodge(width = 0.6)
  ) + 
  
  # prettys 
  scale_color_manual(
    values = c("darkgoldenrod2", "royalblue4"), 
    labels = c("C" = "Control", "W" = "Hydration"), 
    name = "Treatment"
  ) +
  scale_shape_manual(
    values = c(15, 17), 
    labels = c("C" = "Control", "W" = "Hydration"),
    name = "Treatment"
  ) +
  scale_y_continuous(
    limits = c(240, 410), 
    breaks = seq(200, 400, by = 50),
    name = expression('Osmolality (mmol '*kg^-1*')'),
  ) + 
  scale_x_discrete(
    labels = c("DAB" = "At Birth", "PS" = "Post-Shed"),
    name = NULL
  ) +
  
  theme(legend.position = "none") +
  
  # statistical labels
  annotate(geom = "text", x = 1, y = 405, label = "NS", size = 2) +
  annotate(geom = "text", x = 2, y = 400, label = "*", size = 6) +
  annotate(geom = "text", x = 2.5, y = 335, label = "*", size = 6, 
           color = "darkgoldenrod2") +
  annotate(geom = "text", x = 2.5, y = 290, label = "NS", size = 2,
           color = "royalblue4") -> neo_osmol

neo_osmol

cewl raw

ggplot() +
  geom_pointrange(
    data = dab_ps_compare_means,
    aes(
      x = timept,
      y = CEWL_g_m2h_mean,
      ymin = CEWL_g_m2h_mean-CEWL_g_m2h_sd,
      ymax = CEWL_g_m2h_mean+CEWL_g_m2h_sd,
      group = Treatment,
      color = Treatment,
      shape = Treatment
    ),
    position = position_jitterdodge(jitter.width = 0.6, dodge.width = 0.6),
    size = 0.2, linewidth = 0.3,
    alpha = 0.6,
    show.legend = FALSE
  ) +
  geom_line(
    data = cewl_lm_time_tmt_CI,
    aes(
      x = timept, y = estimate, 
      group = Treatment, 
      color = Treatment
    ),
    position = position_dodge(width = 0.6)
  ) +
  geom_errorbar(
    data = cewl_lm_time_tmt_CI,
    aes(
      x = timept, y = estimate, 
      ymin = conf.low, ymax = conf.high,  
      group = Treatment, 
      color = Treatment
    ),
    width = 0.4,
    position = position_dodge(width = 0.6)
  ) +
  geom_point(
    data = cewl_lm_time_tmt_CI,
    aes(
      x = timept, y = estimate, 
      group = Treatment, 
      color = Treatment,
      shape = Treatment
    ),
    size = 2.5,
    position = position_dodge(width = 0.6)
  ) + 
  
  # prettys 
  scale_color_manual(
    values = c("darkgoldenrod2", "royalblue4"), 
    labels = c("C" = "Control", "W" = "Hydration"), 
    name = "Treatment"
  ) +
  scale_shape_manual(
    values = c(15, 17), 
    labels = c("C" = "Control", "W" = "Hydration"),
    name = "Treatment"
  ) +
  scale_y_continuous(
    limits = c(5, 23), 
    breaks = seq(5, 20, by = 5),
    name = expression('CEWL (g '*m^-2*' '*h^-1*')')
  ) + 
  scale_x_discrete(
    labels = c("DAB" = "At Birth", "PS" = "Post-Shed"),
    name = NULL
  ) +
  
  theme(legend.position = "none") +
  
  # statistical labels
  annotate(geom = "text", x = 1, y = 22.5, label = "NS", size = 2) +
  annotate(geom = "text", x = 2, y = 22.5, label = "NS", size = 2) +
  annotate(geom = "text", x = 2.45, y = 15.5, label = "*", size = 6, 
           color = "darkgoldenrod2") +
  annotate(geom = "text", x = 2.55, y = 15.5, label = "*", size = 6,
           color = "royalblue4") -> neo_cewl
  
neo_cewl

change in mass

ggplot(
  data = deltas_befaft, 
  aes(
    x = timept,
    y = delta_mass,
    group = Mother_ID,
    color = Treatment,
    shape = Treatment
  )) +
  geom_hline(yintercept = 0, linetype = 'dashed', color = "darkgray") +
  geom_point(size = 1.4, position = position_dodge(width = 0.1)) +
  geom_line(alpha = 0.4, position = position_dodge(width = 0.1)) +
  
  # prettys 
  scale_color_manual(
    values = c("darkgoldenrod2", "royalblue4"), 
    labels = c("C" = "Control", "W" = "Hydration"), 
    name = NULL
  ) +
  scale_shape_manual(
    values = c(15, 17), 
    labels = c("C" = "Control", "W" = "Hydration"),
    name = NULL
  ) +
  scale_y_continuous(
    breaks = c(-3, -2, -1, 0),
    name = expression(Delta ~ 'Mass (g)')
  ) +
  scale_x_discrete(
    labels = c("DAB" = "At Birth", "PS" = "Post-Shed"),
    name = NULL
  ) +
  
  theme(legend.position = "inside",
        legend.position.inside = c(0.25, 0.15),
        legend.background = element_blank()) +
  
  # labels
  annotate(geom = "text", x = 2.3, y = -1.1, label = "*", size = 6, 
           color = "royalblue4") +
  annotate(geom = "text", x = 2.3, y = -2.1, label = "*", size = 6, 
           color = "darkgoldenrod2") +
  annotate(geom = "text", x = 2, y = -0.5, label = "*", 
           size = 6) -> neo_change_mass
neo_change_mass

change in osmol

ggplot(
  data = deltas_befaft, 
  aes(
    x = timept,
    y = delta_osml,
    group = Mother_ID,
    color = Treatment,
    shape = Treatment
  )) +
  geom_hline(yintercept = 0, linetype = 'dashed', color = "darkgray") +
  geom_point(size = 1.4, position = position_dodge(width = 0.1)) +
  geom_line(alpha = 0.4, position = position_dodge(width = 0.1)) +
  
  # prettys 
  scale_color_manual(
    values = c("darkgoldenrod2", "royalblue4"), 
    labels = c("C" = "Control", "W" = "Hydration"), 
    name = NULL
  ) +
  scale_shape_manual(
    values = c(15, 17), 
    labels = c("C" = "Control", "W" = "Hydration"),
    name = NULL
  ) +
  scale_y_continuous(
    breaks = c(-20, 0, 20, 40, 60, 80),
    name = expression(Delta ~ 'Osmolality (mmol '*kg^-1*')'),
  ) +
  scale_x_discrete(
    labels = c("DAB" = "At Birth", "PS" = "Post-Shed"),
    name = NULL
  ) +
  
  theme(legend.position = "inside",
        legend.position.inside = c(0.25, 0.85),
        legend.background = element_blank()) +
  
  # labels
  annotate(geom = "text", x = 2.3, y = -5, label = "NS", size = 2, 
           color = "royalblue4") +
  annotate(geom = "text", x = 2.3, y = 45, label = "*", size = 6, 
           color = "darkgoldenrod2") +
  annotate(geom = "text", x = 2, y = 81, label = "*", 
           size = 6) -> neo_change_osml
neo_change_osml

change in cewl

ggplot(
  data = deltas_befaft, 
  aes(
    x = timept,
    y = delta_CEWL,
    group = Mother_ID,
    color = Treatment,
    shape = Treatment
  )) +
  geom_hline(yintercept = 0, linetype = 'dashed', color = "darkgray") +
  geom_point(size = 1.4, position = position_dodge(width = 0.1)) +
  geom_line(alpha = 0.4, position = position_dodge(width = 0.1)) +
  
  # prettys 
  scale_color_manual(
    values = c("darkgoldenrod2", "royalblue4"), 
    labels = c("C" = "Control", "W" = "Hydration"), 
    name = NULL
  ) +
  scale_shape_manual(
    values = c(15, 17), 
    labels = c("C" = "Control", "W" = "Hydration"),
    name = NULL
  ) +
  scale_y_continuous(
    breaks = c(0, 5, 10),
    name = expression(Delta ~ 'CEWL (g '*m^-2*' '*h^-1*')')
  ) +
  scale_x_discrete(
    labels = c("DAB" = "At Birth", "PS" = "Post-Shed"),
    name = NULL
  ) +
  
  theme(legend.position = "inside",
        legend.position.inside = c(0.25, 0.85),
        legend.background = element_blank()) +
  
  # labels
  annotate(geom = "text", x = 2.25, y = 5, label = "*", size = 6,
           color = "darkgoldenrod2") +
  annotate(geom = "text", x = 2.35, y = 5, label = "*", size = 6,
           color = "royalblue4") +
  annotate(geom = "text", x = 2, y = 11, label = "NS", 
           size = 2) -> neo_change_cewl
neo_change_cewl

arrange

ggarrange(
  neo_mass, neo_change_mass,
  neo_osmol, neo_change_osml,
  neo_cewl, neo_change_cewl,
  nrow = 3,
  ncol = 2,
  labels = c("A", "B", "C", "D", "E", "F"),
  align = "hv"
) -> neo_tmt_effects

ggsave(
  filename = "Fig6_neonate_tmt_effects.pdf",
  plot = neo_tmt_effects,
  path = "./results_figures",
  device = "pdf",
  dpi = 600,
  units = "mm",
  width = 160,
  height = 180
)

ggsave(
  filename = "Fig6_neonate_tmt_effects.tiff",
  plot = neo_tmt_effects,
  path = "./results_figures",
  device = "tiff",
  dpi = 600,
  units = "mm",
  width = 160,
  height = 180
)